在R语言中带约束优化多元函数S.residuum的实现咨询

Hey there! Let's tackle this constrained optimization problem in R. Since you're working with a 6-dimensional weight vector w and need to add two missing constraints (one involving stdv and corr), I'll walk you through the most common tools and examples tailored to typical scenarios (I'll assume the stdv/corr constraint is a quadratic risk constraint—super common in quantitative tasks, but you can adjust it to your actual needs).

First: Clarify Your Constraint Form #

Before diving into code, you need to translate your constraints into mathematical terms that R's optimizers can understand:

• For linear constraints (like sum of weights equals 1), base R's constrOptim() works great.
• For nonlinear constraints (like variance limits using stdv and corr), the nloptr package is more flexible.

Let's cover both cases with concrete examples.

Method 1: Using constrOptim() for Linear Constraints #

If your missing constraints are linear (e.g., sum of weights equals 1, or a linear combination involving stdv), this is the simplest approach. constrOptim() requires constraints in the form A %*% w >= b.

Example Code #

Suppose your constraints are:

1. Sum of weights equals 1 (we split this into two inequalities: sum(w) >=1 and sum(w) <=1)
2. A linear constraint using stdv: w %*% stdv >= 0.3 (replace this with your actual linear constraint)

# First, define your target function (replace with your actual S.residuum)
S.residuum <- function(w) {
  # Example: Minimize squared difference from a target vector
  sum((w - c(2, 3, 1, 4, 0, 2))^2)
}

# Initial weight vector
w0 <- c(1, 1, 1, 1, 1, 1)

# Define your stdv vector (replace with your actual values)
stdv <- c(0.1, 0.2, 0.15, 0.3, 0.05, 0.25)

# Build constraint matrix A and vector b for A %*% w >= b
A <- rbind(
  rep(1, 6),          # Constraint 1: sum(w) >= 1
  rep(-1, 6),         # Constraint 2: -sum(w) >= -1 → sum(w) <= 1
  stdv                # Constraint 3: w %*% stdv >= 0.3 (your stdv-related constraint)
)
b <- c(1, -1, 0.3)

# Run constrained optimization
optim_result <- constrOptim(
  theta = w0,
  f = S.residuum,
  grad = NULL,  # Let R compute numerical gradient if you don't have an analytical one
  ui = A,
  ci = b
)

# Check results
cat("Optimal weights:", optim_result$par, "\n")
cat("Minimum S.residuum value:", optim_result$value, "\n")

Method 2: Using nloptr() for Nonlinear Constraints #

If your stdv/corr constraint is nonlinear (e.g., limiting the variance of the weight portfolio: t(w) %*% diag(stdv) %*% corr %*% diag(stdv) %*% w <= target_var), nloptr is the way to go—it supports both linear and nonlinear equality/inequality constraints.

Example Code #

Let's assume your missing constraints are:

1. Nonlinear inequality: Portfolio variance ≤ 0.02 (using stdv and corr)
2. Linear equality: Sum of weights = 1
3. Boundary constraint: All weights ≥ 0 (no shorting)

# Install and load nloptr if you haven't already
# install.packages("nloptr")
library(nloptr)

# Define your target function
eval_f <- function(w) {
  return(S.residuum(w))
}

# Define equality constraint: sum(w) - 1 = 0
eval_g_eq <- function(w) {
  return(sum(w) - 1)
}

# Define your stdv and corr matrix (replace with your actual values)
stdv <- c(0.1, 0.2, 0.15, 0.3, 0.05, 0.25)
corr <- matrix(runif(6*6), nrow=6)
corr <- (corr + t(corr))/2  # Ensure correlation matrix is symmetric
diag(corr) <- 1             # Diagonal elements are 1
cov_mat <- diag(stdv) %*% corr %*% diag(stdv)  # Covariance matrix
target_var <- 0.02          # Target maximum variance

# Define nonlinear inequality constraint: variance - target_var ≤ 0
eval_g_ineq <- function(w) {
  return(t(w) %*% cov_mat %*% w - target_var)
}

# Boundary constraints: weights ≥ 0, no upper limit (adjust ub if needed)
lb <- rep(0, 6)
ub <- rep(Inf, 6)

# Initial weights
w0 <- c(1, 1, 1, 1, 1, 1)

# Optimization options (SLSQP algorithm works well for mixed constraints)
opts <- list(
  algorithm = "NLOPT_LD_SLSQP",
  xtol_rel = 1e-6,
  maxeval = 1000
)

# Run optimization
optim_result <- nloptr(
  x0 = w0,
  eval_f = eval_f,
  eval_g_eq = eval_g_eq,
  eval_g_ineq = eval_g_ineq,
  lb = lb,
  ub = ub,
  opts = opts
)

# Check results
cat("Optimal weights:", optim_result$solution, "\n")
cat("Minimum S.residuum value:", optim_result$objective, "\n")

Quick Tips for Success #

• If you can compute the analytical gradient of S.residuum, pass it to the optimizer (e.g., grad in constrOptim() or eval_grad_f in nloptr)—it'll make the optimization faster and more reliable.
• Double-check constraint directions: constrOptim() expects ui %*% w >= ci, while nloptr's inequality constraints are eval_g_ineq(w) <= 0.
• Make sure your initial w0 is in the feasible region (satisfies all constraints)—otherwise the optimizer might struggle to converge.

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