optimization - 使用 RcppEnsmallen 发散的分位数回归

Question

我正在尝试RcppEnsmallen 估计分位数回归模型。我注意到结果似乎并不收敛。每次弹出不同的结果。我在我的代码中设置了它的一阶条件。我想这种分歧可能是由于目标函数在零处不可微。我不确定。

这里的代码

#include <RcppEnsmallen.h>
// [[Rcpp::depends(RcppEnsmallen)]]

class QuantileRegressionFunction
{
public:
  // Construct the QuantileRegressFunction with the given data.
  QuantileRegressionFunction(const arma::mat& X,
                             const arma::colvec& y,
                             const double tau
                  ) : X(X),y(y),tau(tau){}
  
  // Define the objective function.
  double EvaluateWithGradient(const arma::mat& beta, arma::mat& gradient ){
    const arma::mat fit = y - X * beta;
    const arma::vec eval =0.5* arma::abs(fit) + (tau - 0.5  )*fit;
    const arma::vec v= tau - arma::ones<arma::vec>(y.n_rows) % ( fit < 0 )  ;
    gradient = - (X.t() * v)*(1/X.n_rows) ;
    return   arma::accu( (1/y.n_rows) *eval ) ;
    
    }
  
  
private:
   const arma::mat& X;
   const arma::vec& y;
   const double tau;
  };

// [[Rcpp::export]]

arma::mat qr_reg(const arma::mat& X, const arma::vec& y,const double& tau){
  QuantileRegressionFunction qrf(X,y,tau);
 
  
  ens::L_BFGS lbfgs;
  lbfgs.MaxIterations()= 1000 ;
  //lbfgs.MaxLineSearchTrials() = 10000;
  //lbfgs.ArmijoConstant() = 1e-18;
  arma::mat beta(X.n_cols, 1 , arma::fill::randn);
  
  lbfgs.Optimize(qrf,beta);
  return beta;
  }

这是一个简单的模拟练习：

 n <- 1000
  beta <- c(-2, 1.5, 3, 8.2, 6.6)
  p <- length(beta)
  X <- cbind(1, matrix(rnorm(n), ncol = p - 1))
  y <- X %*% beta + rnorm(n / (p - 1))
qr_reg(X,y,0.1)

每次qr_reg(X,y,.1)报告不同的估计。

标签： optimizationrcppconvex-optimization