带约束的矩阵线性组合系数最小二乘优化的Python高效实现

你需要找到满足约束 (x_1 + x_2 + x_3 = 1) 且 (0 \leq x_1,x_2,x_3 \leq 1) 的系数 (x_1,x_2,x_3)，最小化目标函数 (|D - x_1A - x_2B - x_3C|_2^2)（即所有元素误差的平方和）。由于变量维度极低（仅2个独立变量，因为 (x_3=1-x_1-x_2)），可以用解析+边界检查的高效方法，或者直接调用优化库求解。

方法一：解析解法（高效，适合低维场景） #

核心思路 #

通过变量替换将问题转化为二元无约束最小二乘，再检查解是否符合约束；若不符合，则在可行域边界上求解一维最小二乘，最终选择损失最小的解。

步骤与代码实现 #

import numpy as np
from scipy.linalg import solve

def find_optimal_coefficients(A, B, C, D):
    # 展平高维矩阵为一维向量
    a = A.flatten()
    b = B.flatten()
    c = C.flatten()
    d = D.flatten()
    
    # 变量替换简化目标函数
    y = d - c
    u = a - c
    v = b - c
    
    # 求解无约束最小二乘的线性方程组
    mat = np.array([
        [np.dot(u, u), np.dot(u, v)],
        [np.dot(v, u), np.dot(v, v)]
    ])
    rhs = np.array([np.dot(u, y), np.dot(v, y)])
    
    # 处理矩阵奇异的特殊情况
    try:
        x1_uc, x2_uc = solve(mat, rhs)
    except np.linalg.LinAlgError:
        x1_uc = 0.0
        x2_uc = np.dot(v, y) / np.dot(v, v) if np.dot(v, v) != 0 else 0.0
    
    x3_uc = 1 - x1_uc - x2_uc
    
    # 检查无约束解是否在可行域内
    if x1_uc >= -1e-8 and x2_uc >= -1e-8 and x3_uc >= -1e-8:
        return (x1_uc, x2_uc, x3_uc)
    
    # 定义损失计算函数
    def compute_loss(x1, x2):
        residual = y - x1*u - x2*v
        return np.dot(residual, residual)
    
    # 边界1：x1=0，求解最优x2
    x2_b1 = np.dot(v, y) / np.dot(v, v) if np.dot(v, v) != 0 else 0.0
    x2_b1 = np.clip(x2_b1, 0, 1)
    loss_b1 = compute_loss(0.0, x2_b1)
    
    # 边界2：x2=0，求解最优x1
    x1_b2 = np.dot(u, y) / np.dot(u, u) if np.dot(u, u) != 0 else 0.0
    x1_b2 = np.clip(x1_b2, 0, 1)
    loss_b2 = compute_loss(x1_b2, 0.0)
    
    # 边界3：x1+x2=1，求解最优t=x1
    w = u - v
    z = y - v
    t_b3 = np.dot(w, z) / np.dot(w, w) if np.dot(w, w) != 0 else 0.0
    t_b3 = np.clip(t_b3, 0, 1)
    loss_b3 = compute_loss(t_b3, 1 - t_b3)
    
    # 选择损失最小的边界解
    min_loss_idx = np.argmin([loss_b1, loss_b2, loss_b3])
    if min_loss_idx == 0:
        return (0.0, x2_b1, 1 - x2_b1)
    elif min_loss_idx == 1:
        return (x1_b2, 0.0, 1 - x1_b2)
    else:
        return (t_b3, 1 - t_b3, 0.0)

# 测试示例
if __name__ == "__main__":
    m, n, p = 2, 2, 2
    A = np.random.rand(m, n, p)
    B = np.random.rand(m, n, p)
    C = np.random.rand(m, n, p)
    D = 0.3*A + 0.5*B + 0.2*C + np.random.randn(m, n, p)*0.01  # 带噪声的目标矩阵
    
    x1, x2, x3 = find_optimal_coefficients(A, B, C, D)
    print(f"最优系数: x1={x1:.4f}, x2={x2:.4f}, x3={x3:.4f}")
    print(f"约束验证: x1+x2+x3={x1+x2+x3:.4f}, 系数非负: {x1>=0 and x2>=0 and x3>=0}")

方法二：调用优化库（简洁，无需手动推导） #

如果不想处理复杂的边界逻辑，可以直接使用scipy.optimize.minimize，选择支持约束和边界的SLSQP优化器。

代码实现 #

import numpy as np
from scipy.optimize import minimize

def find_optimal_coefficients_opt(A, B, C, D):
    a = A.flatten()
    b = B.flatten()
    c = C.flatten()
    d = D.flatten()
    
    # 定义目标函数
    def objective(x):
        x1, x2 = x
        x3 = 1 - x1 - x2
        residual = d - x1*a - x2*b - x3*c
        return np.dot(residual, residual)
    
    # 设置约束与边界
    constraints = [{'type': 'ineq', 'fun': lambda x: 1 - x[0] - x[1]}]
    bounds = [(0, 1), (0, 1)]
    
    # 初始猜测值
    x0 = np.array([1/3, 1/3])
    
    # 求解优化问题
    result = minimize(objective, x0, method='SLSQP', bounds=bounds, constraints=constraints)
    
    x1_opt, x2_opt = result.x
    x3_opt = 1 - x1_opt - x2_opt
    return (x1_opt, x2_opt, x3_opt)

# 测试示例
if __name__ == "__main__":
    m, n, p = 2, 2, 2
    A = np.random.rand(m, n, p)
    B = np.random.rand(m, n, p)
    C = np.random.rand(m, n, p)
    D = 0.3*A + 0.5*B + 0.2*C + np.random.randn(m, n, p)*0.01
    
    x1, x2, x3 = find_optimal_coefficients_opt(A, B, C, D)
    print(f"最优系数: x1={x1:.4f}, x2={x2:.4f}, x3={x3:.4f}")
    print(f"约束验证: x1+x2+x3={x1+x2+x3:.4f}, 系数非负: {x1>=0 and x2>=0 and x3>=0}")

方法对比 #

• 解析法：速度极快，仅涉及向量点积和解2x2方程组，无迭代过程，适合大规模矩阵。
• 优化库法：代码简洁，无需手动推导边界逻辑，适合快速实现，对于常规规模矩阵，速度差异可忽略。

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