如何解决REBayes中KWDual的MSK_RES_TRM_STALL优化停滞错误？

Hey, let's tackle this Mosek stall error you're hitting with your R simulation. The MSK_RES_TRM_STALL error usually pops up when the solver can't make meaningful progress, often due to poor scaling, ill-conditioned problems, or near-feasible/infeasible models. Since reducing the simulation size isn't ideal, here are targeted fixes tailored to your code:

1. Fix Model Scaling (Most Impactful) #

Your simulation uses large datasets (n=1000) which creates a huge matrix A (1000x300) in the BDE function, leading to extreme numerical ranges that confuse Mosek. Standardizing your input data will normalize these ranges and help the solver converge.

Modify your sim1 function to scale the data before processing, then reverse the scaling for your final estimates:

sim1 <- function(n, R = 10, setting = 0){
  A <- matrix(0, 4, R)
  if(setting == 0){
    G0 <- function(t) punif(t,0,6)/8 + 7 * pnorm(t, 0, 0.5)/8
    rf0 <- function(n){
      s <- sample(0:1, n, replace = TRUE, prob = c(1,7)/8)
      rnorm(n) + (1-s) * runif(n,0,6) + s * rnorm(n,0,0.5)
    }
  } else{
    G0 <- function(t) 0 + 7 * (t > 0)/8 + (t > 2)/8
    rf0 <- function(n){
      s <- sample(0:1, n, replace = TRUE, prob = c(1,7)/8)
      rnorm(n) + (1-s) * 2 + s * 0
    }
  }
  for(i in 1:R){
    y <- rf0(n)
    # Add data scaling here
    y_mean <- mean(y)
    y_sd <- sd(y)
    y_scaled <- (y - y_mean)/y_sd
    
    # Use scaled y for all estimators
    g <- BDE(y_scaled)
    # Reverse scaling for density estimates
    g$x <- g$x * y_sd + y_mean
    g$y <- g$y / y_sd  # Density scaling rule: divide by sd
    Wg <- Wasser(G0, g)
    
    h <- GLmix(y_scaled)
    h$x <- h$x * y_sd + y_mean
    h$y <- h$y / y_sd
    Wh <- Wasser(G0, h)
    Whs <- Wasser(G0, h, interp = "biweight")
    
    k <- KFE(y_scaled)
    k$x <- k$x * y_sd + y_mean
    k$y <- k$y / y_sd
    Wk <- Wasser(G0, k)
    
    A[,i] <- c(Wg$W, Wk$W, Wh$W, Whs$W)
  }
  A
}

2. Tune Mosek Solver Parameters #

You can pass custom Mosek parameters through the REBayes package's GLmix function to adjust the solver's behavior and avoid stalls. Key parameters to tweak include iteration limits, stall tolerance, and scaling mode:

# Modify the GLmix call in sim1 to include Mosek parameters
h <- GLmix(y_scaled, control = list(
  mskparam = list(
    MSK_IPAR_INTPNT_MAX_ITERATIONS = 1000,  # Increase max iterations
    MSK_IPAR_INTPNT_TOL_STALL = 1e-6,       # Lower stall tolerance for stricter progress checks
    MSK_IPAR_INTPNT_SCALING = MSK_SCALING_MODE_AGGRESSIVE  # Enable aggressive scaling
  )
))

Note: These Mosek parameter constants are included with the RMosek package, so ensure it's properly installed and loaded.

3. Simplify Model Complexity #

Reducing the complexity of your estimation models can make the problem more numerically stable:

• In the BDE function, lower the spline degrees of freedom (e.g., df=3 instead of df=5) to reduce the number of variables in the optimization:

  BDE <- function(y, T = 300, df = 3, c0 = 1){  # Changed df from 5 to 3
    # ... rest of the function remains the same
  }
  

• For GLmix, reduce the number of mixture components with the K argument (the default may be too high for large datasets):

  h <- GLmix(y_scaled, K=20)  # Lower K to reduce model complexity
  

4. Fix Numerical Overflow in BDE #

The exp(X %*% a) in your qmle function can cause numerical overflow with large datasets, leading to ill-conditioned objective functions. Add a log-sum-exp-style adjustment to stabilize calculations:

qmle <- function(a, X, A, c0){
  linear_pred <- X %*% a
  linear_pred <- linear_pred - max(linear_pred)  # Prevent exp overflow by centering
  g <- exp(linear_pred)
  g <- g/sum(g)
  f <- A %*% g
  -sum(log(f)) + c0 * sum(a^2)^.5
}

Start with the scaling fix first—it's the most common solution for MSK_RES_TRM_STALL errors. If that doesn't work, try combining it with the parameter tuning and model simplification steps.

内容的提问来源于stack exchange，提问作者Alex