如何在Ax >= b约束下求解带非负x的最小二乘解？

阿华AIGC实验室

2026-4-15

嗨，我明白你的问题了——你用scipy.optimize.nnls得到的结果里有些营养没达标（Ax的元素小于b），这是因为nnls的核心目标是让Ax尽可能接近b（最小化平方误差），但它不会强制Ax的每一项都大于等于b，甚至为了让整体误差更小，可能会牺牲某些维度的达标情况。

你的需求本质上是一个带约束的凸优化问题：要找非负的x，让Ax的每一项都不小于b（满足所有营养摄入需求），同时最好能让某个目标最优（比如总食物量最少，或者多余的营养最少）。这种问题用普通的最小二乘工具搞不定，得用**线性规划（LP）**或者专门的凸优化库来解决，下面给你两种具体的实现方法：

方法一：用cvxpy（直观易写，推荐）

cvxpy是专门用于凸优化的Python库，它允许你用接近数学表达式的方式定义目标函数和约束，非常适合你的营养问题。比如你可以轻松定义“总食物量最少”或者“多余营养最少”的目标，同时加上Ax≥b和x≥0的约束。

代码示例

import cvxpy as cp
import numpy as np

# 复用你提供的营养矩阵和需求向量
v1 = np.array([0.00000961,0.0000011,0.0000011,0.000015,0.00000884,0.00000286,0,0.00000006,0,0.000196,0,0.0000071,0.000000023,0.000131,0.00038,0,0,0.00000161,0,0,0.0000069,0.00027,0.000005,0,0,0.00054,0.00475,0.000000002,0.00036,0.0000032,0,0,0.033,0.0015,0.02,0.00052,0.207]).T
v2 = np.array([0,0.0000064,0.00000135,0.000121,0.0000177,0.00000348,0,0.00000006,0,0,0,0.0000833,0,0.000525,0.00062,0,0,0.0000114,0,0,0.0000458,0.00168,0.0000193,0,0,0.00376,0.00705,0.000000072,0.00018,0.0000327,0,0,0.085,0.492,0.258,0.0628,0.161]).T
v3 = np.array([0,0.000009,0.00000193,0.0000196,0.00000899,0,0,0.00000444,0,0,0,0.0000021,0.000000056,0.000664,0.00123,0,0,0.00000841,0,0,0.0000502,0.00171,0.0000106,0,0,0.00352,0.0148,0.000000032,0.00005,0.0000365,0,0,0.155,0.0142,0.216,0.00366,0.624]).T
v4 = np.array([0,0.00000763,0.00000139,0.00000961,0.0000135,0.00000119,0,0.00000056,0,0,0,0,0,0,0.00054,0,0,0.00000626,0,0,0.0000472,0.00177,0.0000492,0,0,0.00523,0.00429,0,0.00002,0.0000397,0,0,0.106,0.069,0.169,0.0122,0.663]).T
v5 = np.array([0.00000241,0.00000113,0.00000347,0.0000118,0.0000037,0.00000147,0,0.00000062,0,0.000934,0,0,0,0.000005,0.00254,0,0,0.00000053,0,0,0.000016,0.00033,0.0000092,0,0,0.00055,0.00348,0,0.00053,0.0000039,0,0,0.041,0.0149,0.0292,0.00178,0.0442]).T
b = np.array([0.00657,0.00876,0.00949,0.1168,0.0365,0.01241,0.000219,0.00292,0,0.657,0.000146,0.1095,0.000876,4.015,8.76,16.79,0.0002555,0.00657,0.0292,0.001095,0.0584,3.066,0.01679,0.0003285,0,5.11,24.82,0.0004015,10.95,0.0803,0,0.00219,204.4,569.4,365,146,2007.5]).T
A = np.column_stack((v1,v2,v3,v4,v5))

# 定义优化变量x：5种食物的量，强制非负
x = cp.Variable(5, nonneg=True)

# 约束条件：每种营养摄入都达标（Ax >= b）
constraints = [A @ x >= b]

# 选择目标函数：这里选「总食物量最少」，符合实际饮食需求
# 若想追求营养刚好达标，可替换为 cp.Minimize(cp.sum(A @ x - b))
objective = cp.Minimize(cp.sum(x))

# 构建问题并求解
prob = cp.Problem(objective, constraints)
prob.solve(solver=cp.ECOS)

# 输出结果
print("最优食物量x：")
print(x.value)
print("\n实际摄入营养Ax：")
print(A @ x.value)
print("\n营养达标情况（Ax - b，所有值应≥0）：")
print((A @ x.value) - b)

方法二：用scipy.optimize.linprog（无需额外安装库）

如果你不想安装新库，也可以用scipy自带的linprog工具，它专门解决线性规划问题。需要注意linprog默认支持A_ub @ x <= b_ub形式的约束，因此要把Ax >= b转化为-Ax <= -b。

代码示例

import numpy as np
from scipy.optimize import linprog

# 复用你提供的营养矩阵和需求向量（和上面一致，此处省略重复代码）
# ...

# 目标函数：最小化总食物量，系数为全1向量
c = np.ones(5)

# 约束转换：Ax >= b → -Ax <= -b
A_ub = -A
b_ub = -b

# 变量边界：每种食物的量不能为负
bounds = [(0, None) for _ in range(5)]

# 求解线性规划，推荐使用highs求解器（scipy 1.9+版本支持）
res = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=bounds, method='highs')

# 输出结果
print("最优食物量x：")
print(res.x)
print("\n实际摄入营养Ax：")
print(A @ res.x)
print("\n营养达标情况（Ax - b，所有值应≥0）：")
print((A @ res.x) - b)

关键说明

为什么nnls不适用？
scipy.optimize.nnls的核心目标是最小化||Ax - b||₂²（Ax与b的平方误差），同时保证x≥0。但这个目标允许Ax在某些维度上小于b——因为减小这些维度的Ax可能让整体平方误差更小。而你的需求是所有营养维度必须达标，所以nnls无法满足约束。
约束可行性检查
如果求解器返回“infeasible”，说明当前的5种食物无法组合出满足需求b的营养摄入，你可能需要调整营养需求b，或者增加新的食物种类。
目标函数的灵活选择
- 最小化总食物量：适合追求“少吃点”的场景
- 最小化多余营养：适合追求“刚好达标、不浪费”的场景
- 最小化x的L2范数：适合让食物种类的摄入量更均衡的场景

备注：内容来源于stack exchange，提问作者Arcaeca