带约束的最小二乘法
least-squares method with a constraint
我有 37 个具有 36 个变量的线性方程，其形式为矩阵：A x = b。 （A 有 37 行和 36 列。）方程没有精确解，所以我使用 Matlab 使用 x = A \ b.
问题是我还有一个条件，x的所有元素都应该是正数：xi > 0 for all i。 x = A \ b 对某些元素给出负值。我怎样才能应用这个约束？
以下是我正在使用的 A 和 b 的具体值：
A = [0.83 0.17 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0.02 0.63 0.17 0 0 0 0.18 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0.02 -0.37 0.17 0 0 0 0.18 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0.02 -0.37 0.17 0 0 0 0.18 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0.02 -0.37 0.17 0 0 0 0.18 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0.02 -0.2 0 0 0 0 0.18 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0.15 0 0 0 0 0 -0.32 0.17 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0.15 0 0 0 0 0.02 -0.52 0.17 0 0 0 0.18 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0.15 0 0 0 0 0.02 -0.52 0.17 0 0 0 0.18 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0.15 0 0 0 0 0.02 -0.52 0.17 0 0 0 0.18 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0.15 0 0 0 0 0.02 -0.52 0.17 0 0 0 0.18 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0.15 0 0 0 0 0.02 -0.33 0 0 0 0 0.18 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0.15 0 0 0 0 0 -0.32 0.17 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0.15 0 0 0 0 0.02 -0.52 0.17 0 0 0 0.18 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0.15 0 0 0 0 0.02 -0.52 0.17 0 0 0 0.18 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0.02 -0.52 0.17 0 0 0 0.18 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0.02 -0.52 0.17 0 0 0 0.18 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0.02 -0.35 0 0 0 0 0.18 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0 -0.32 0.17 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0.02 -0.52 0.17 0 0 0 0.18 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0.02 -0.52 0.17 0 0 0 0.18 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0.02 -0.52 0.17 0 0 0 0.18 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0.02 -0.52 0.17 0 0 0 0.18 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0.02 -0.35 0 0 0 0 0.18 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0 -0.32 0.17 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0.02 -0.52 0.17 0 0 0 0.18 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0.02 -0.52 0.17 0 0 0 0.18 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0.02 -0.52 0.17 0 0 0 0.18 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0.02 -0.52 0.17 0 0 0 0.18 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0.02 -0.35 0 0 0 0 0.018 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0 -0.32 0.17 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0.02 -0.34 0.17 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0.02 -0.34 0.17 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0.02 -0.34 0.17 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0.02 -0.34 0.17
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0.02 -0.17
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1];
b = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1]';
当 x>0 条件时，这将是一个线性优化问题。解决该问题的最佳算法是单纯形算法。这个想法是每个线性方程 aix=bi 提供一条线，这些线的组合提供一个 polygon/polyhedron，答案是这个 polyhedron/polygon 的顶点之一。单纯形算法非常标准，有许多可用的函数和库可以计算它。
您想在最小二乘意义上找到 A x = b 的近似解 x，即您想要最小化
||A x - b||2 = xT AT A x
+ bT b - 2 xT AT b.
不计常数项bTb除以2，符合aquadratic programming problem的形式，即最小化
1/2 xT H x + fT x,
如果我们选择H = AT A and f = - AT b.
quadprog对应的用法是：
H = A' * A;
f = - A' * b;
x = quadprog(H, f)
您还希望 x 的元素为正。可以使用 quadprog、A 和 b 的附加参数引入非负约束（不要与您的矩阵混淆！）：
n = size(A, 2);
x = quadprog(H, f, -eye(n), zeros(n, 1))
正性约束没有意义，因为如果最优解涉及x的一个或多个元素恰好为0，那么严格正解会更好对应元素越小：0.01会更好比 0.1，0.001 会比 0.01 好，等等——没有自然界。如果你想确保解决方案是全正的，你必须自己设置一个有限界限：
x = quadprog(H, f, -eye(n), zeros(n, 1) + 0.001)
现在 x 的元素的最小可能值为 0.001。
更新问题后补充了A和b的实际数据：使用代码
H = A' * A;
f = - A' * b;
n = size(A, 2);
x = quadprog(H, f, -eye(n), zeros(n, 1))
我得到结果：
Minimum found that satisfies the constraints.
Optimization completed because the objective function is non-decreasing in
feasible directions, to within the default value of the function tolerance,
and constraints are satisfied to within the default value of the constraint tolerance.
<stopping criteria details>
x =
0.000380906335150292
3.90638261088393e-05
0.0111196970167585
0.0227055107206744
0.0318402514628274
0.0371743514880516
0.000800900221354844
0.00746652476710186
0.0180511534370576
0.0282423767946842
0.0362606972021829
0.0417582260990786
0.00860220929402253
0.0174105435824309
0.0265771677458008
0.0343071472371469
0.0395176470725881
0.0419494410289298
0.0187719294637544
0.0268976053211278
0.0336818044612046
0.0382365751296441
0.0398823076542831
0.0391016682549663
0.0279383031707377
0.0339393563379992
0.0377917413001034
0.0382731422972829
0.0338557405807941
0.0217568643500703
0.0343698083354502
0.0381554349806972
0.0392353941260779
0.0368010570888738
0.031271868401718
0.0258232230013864
Matlab 函数lsqlin 或lsqnonneg 可以用来解决您的问题。
例如：
x=lsqnonneg(A,b)
会给你想要的。
我有 37 个具有 36 个变量的线性方程，其形式为矩阵：A x = b。 （A 有 37 行和 36 列。）方程没有精确解，所以我使用 Matlab x = A \ b.
问题是我还有一个条件，x的所有元素都应该是正数：xi > 0 for all i。 x = A \ b 对某些元素给出负值。我怎样才能应用这个约束？
以下是我正在使用的 A 和 b 的具体值：
A = [0.83 0.17 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0.02 0.63 0.17 0 0 0 0.18 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0.02 -0.37 0.17 0 0 0 0.18 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0.02 -0.37 0.17 0 0 0 0.18 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0.02 -0.37 0.17 0 0 0 0.18 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0.02 -0.2 0 0 0 0 0.18 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0.15 0 0 0 0 0 -0.32 0.17 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0.15 0 0 0 0 0.02 -0.52 0.17 0 0 0 0.18 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0.15 0 0 0 0 0.02 -0.52 0.17 0 0 0 0.18 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0.15 0 0 0 0 0.02 -0.52 0.17 0 0 0 0.18 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0.15 0 0 0 0 0.02 -0.52 0.17 0 0 0 0.18 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0.15 0 0 0 0 0.02 -0.33 0 0 0 0 0.18 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0.15 0 0 0 0 0 -0.32 0.17 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0.15 0 0 0 0 0.02 -0.52 0.17 0 0 0 0.18 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0.15 0 0 0 0 0.02 -0.52 0.17 0 0 0 0.18 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0.02 -0.52 0.17 0 0 0 0.18 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0.02 -0.52 0.17 0 0 0 0.18 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0.02 -0.35 0 0 0 0 0.18 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0 -0.32 0.17 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0.02 -0.52 0.17 0 0 0 0.18 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0.02 -0.52 0.17 0 0 0 0.18 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0.02 -0.52 0.17 0 0 0 0.18 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0.02 -0.52 0.17 0 0 0 0.18 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0.02 -0.35 0 0 0 0 0.18 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0 -0.32 0.17 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0.02 -0.52 0.17 0 0 0 0.18 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0.02 -0.52 0.17 0 0 0 0.18 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0.02 -0.52 0.17 0 0 0 0.18 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0.02 -0.52 0.17 0 0 0 0.18 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0.02 -0.35 0 0 0 0 0.018 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0 -0.32 0.17 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0.02 -0.34 0.17 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0.02 -0.34 0.17 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0.02 -0.34 0.17 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0.02 -0.34 0.17
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.15 0 0 0 0 0.02 -0.17
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1];
b = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1]';
当 x>0 条件时，这将是一个线性优化问题。解决该问题的最佳算法是单纯形算法。这个想法是每个线性方程 aix=bi 提供一条线，这些线的组合提供一个 polygon/polyhedron，答案是这个 polyhedron/polygon 的顶点之一。单纯形算法非常标准，有许多可用的函数和库可以计算它。
您想在最小二乘意义上找到 A x = b 的近似解 x，即您想要最小化
||A x - b||2 = xT AT A x + bT b - 2 xT AT b.
不计常数项bTb除以2，符合aquadratic programming problem的形式，即最小化
1/2 xT H x + fT x,
如果我们选择H = AT A and f = - AT b.
quadprog对应的用法是：
H = A' * A;
f = - A' * b;
x = quadprog(H, f)
您还希望 x 的元素为正。可以使用 quadprog、A 和 b 的附加参数引入非负约束（不要与您的矩阵混淆！）：
n = size(A, 2);
x = quadprog(H, f, -eye(n), zeros(n, 1))
正性约束没有意义，因为如果最优解涉及x的一个或多个元素恰好为0，那么严格正解会更好对应元素越小：0.01会更好比 0.1，0.001 会比 0.01 好，等等——没有自然界。如果你想确保解决方案是全正的，你必须自己设置一个有限界限：
x = quadprog(H, f, -eye(n), zeros(n, 1) + 0.001)
现在 x 的元素的最小可能值为 0.001。
更新问题后补充了A和b的实际数据：使用代码
H = A' * A;
f = - A' * b;
n = size(A, 2);
x = quadprog(H, f, -eye(n), zeros(n, 1))
我得到结果：
Minimum found that satisfies the constraints.
Optimization completed because the objective function is non-decreasing in
feasible directions, to within the default value of the function tolerance,
and constraints are satisfied to within the default value of the constraint tolerance.
<stopping criteria details>
x =
0.000380906335150292
3.90638261088393e-05
0.0111196970167585
0.0227055107206744
0.0318402514628274
0.0371743514880516
0.000800900221354844
0.00746652476710186
0.0180511534370576
0.0282423767946842
0.0362606972021829
0.0417582260990786
0.00860220929402253
0.0174105435824309
0.0265771677458008
0.0343071472371469
0.0395176470725881
0.0419494410289298
0.0187719294637544
0.0268976053211278
0.0336818044612046
0.0382365751296441
0.0398823076542831
0.0391016682549663
0.0279383031707377
0.0339393563379992
0.0377917413001034
0.0382731422972829
0.0338557405807941
0.0217568643500703
0.0343698083354502
0.0381554349806972
0.0392353941260779
0.0368010570888738
0.031271868401718
0.0258232230013864
Matlab 函数lsqlin 或lsqnonneg 可以用来解决您的问题。
例如：
x=lsqnonneg(A,b)
会给你想要的。