如何在R中求解带最小化约束的线性方程组以确定权重？

这是个典型的带线性约束的二次最小化问题，目标是最小化权重的平方和，同时满足两个线性等式约束。我给你分两种方法来解决：一种是通过拉格朗日乘数法推导解析解，另一种是直接用R的优化包做数值求解。

首先明确咱们的问题：

• 给定输入序列 Inputs <- seq(2,7.7,0.3)（共20个元素）
• 要求解权重 w1~w20，满足：
  1. sum(w * Inputs) = 4.8
  2. sum(w) = 1
  3. 最小化目标函数 sum(w^2)

方法一：解析解法（拉格朗日乘数法） #

因为目标函数是凸的，约束是线性的，我们可以用拉格朗日乘数法推导出精确解，不需要迭代计算。

推导过程 #

构造拉格朗日函数：

L = sum(w_i²) - λ(sum(w_i) - 1) - μ(sum(w_i*x_i) - 4.8)

对每个w_i求偏导并令其为0：

∂L/∂w_i = 2w_i - λ - μx_i = 0 → w_i = (λ + μx_i)/2

把这个表达式代入两个约束条件，得到关于λ和μ的二元一次方程组：

1. 代入sum(w_i)=1：
   (20λ + μsum(x_i))/2 = 1 → 20λ + μsum_x = 2
2. 代入sum(w_i*x_i)=4.8：
   (λsum(x_i) + μsum(x_i²))/2 = 4.8 → λsum_x + μsum_x2 = 9.6

解这个方程组就能得到λ和μ，再代入w_i的表达式就能算出所有权重。

R代码实现解析解 #

# 定义输入序列
Inputs <- seq(2,7.7,0.3)
n <- length(Inputs)

# 计算所需的求和项
sum_x <- sum(Inputs)
sum_x2 <- sum(Inputs^2)

# 构造方程组：A %*% c(λ, μ) = b
A <- matrix(c(n, sum_x, sum_x, sum_x2), nrow=2)
b <- c(2, 9.6)

# 解方程组得到λ和μ
lambda_mu <- solve(A, b)
lambda <- lambda_mu[1]
mu <- lambda_mu[2]

# 计算权重
Weights <- (lambda + mu * Inputs)/2

# 验证约束条件
cat("sum(Weights) =", sum(Weights), "\n")
cat("sum(Weights*Inputs) =", sum(Weights*Inputs), "\n")
cat("sum(Weights^2) =", sum(Weights^2), "\n")

方法二：数值解法（用quadprog包） #

如果不想手动推导，直接用R的quadprog包来求解二次规划问题会更快捷，这个包专门处理带线性约束的二次最小化问题。

二次规划的标准形式 #

quadprog的核心函数solve.QP要求问题写成如下形式：

min 0.5 * t(w) %% D %% w + t(d) %% w
约束条件：t(A_eq) %% w = b_eq

对应咱们的问题：

• 目标函数sum(w^2)等价于0.5 * t(w) %*% (2*diag(n)) %*% w，所以D = 2*diag(n)，d = rep(0, n)
• 等式约束：sum(w)=1和sum(w*Inputs)=4.8，所以A_eq = rbind(rep(1, n), Inputs)，b_eq = c(1, 4.8)

R代码实现数值解 #

# 安装并加载quadprog包（如果没装过的话）
# install.packages("quadprog")
library(quadprog)

Inputs <- seq(2,7.7,0.3)
n <- length(Inputs)

# 构造二次规划所需的矩阵和向量
D <- 2 * diag(n)  # 对应目标函数的二次项
d <- rep(0, n)    # 线性项系数为0

# 等式约束
A_eq <- rbind(rep(1, n), Inputs)
b_eq <- c(1, 4.8)

# 求解二次规划
result <- solve.QP(Dmat = D, dvec = d, Amat = t(A_eq), bvec = b_eq, meq = 2)
# meq=2表示前2个约束是等式约束

# 提取权重解
Weights_qp <- result$solution

# 验证结果
cat("sum(Weights_qp) =", sum(Weights_qp), "\n")
cat("sum(Weights_qp*Inputs) =", sum(Weights_qp*Inputs), "\n")
cat("sum(Weights_qp^2) =", sum(Weights_qp^2), "\n")

两种方法得到的结果应该几乎完全一致（数值解法可能有极小的浮点误差），都能满足约束条件并最小化权重的平方和。

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