带约束的最小二乘法答案

【问题标题】：least-squares method with a constraint带约束的最小二乘法
【发布时间】：2015-08-01 23:52:18
【问题描述】：

我有 37 个线性方程，有 36 个矩阵形式的变量：A x = b。（A 有 37 行和 36 列。）方程没有精确解，所以我使用 Matlab to find the closest answer 和 x = A \ b。

问题是我也有一个条件，x的所有元素都应该是正数：x_i > 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]';

【问题讨论】：

如果你有 37 个方程对应 36 个未知数，那是一个超定系统，这意味着不可能有多个解（实际上，很可能没有解）。
我知道，但我不能至少寻找一个错误最少的近似解决方案吗？
如果你想找到最小化错误的X0，首先你需要定义那个错误。 :-) 对于错误E = B - A*X0，您会使用什么类型的规范？欧几里得，出租车，最大值？
@CST-Link，标题说的是最小二乘，这里指定了代价函数
@nasim，谢谢，但这没有意义。这个A是36x36，这个b是1x38。请解决这个问题。

标签： matlab math curve-fitting linear-programming least-squares

【解决方案1】：

Matlab 函数lsqlin 或lsqnonneg 可用于解决您的问题。例如：

x=lsqnonneg(A,b)

会给你你想要的。

【讨论】：

【解决方案2】：

您想在最小二乘意义上找到 A x = b 的近似解 x，即您想最小化

||A x - b||² = x^T A^T A x + b^T b - 2 x^T A^T b.

忽略常数项 b^Tb 并除以因子 2，这符合 quadratic programming problem 的形式，即最小化

1/2 x^T H x + f^T x,

如果我们选择 H = A^T A 和 f = - A^T 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

【讨论】：

警告：Trust-region-reflective 算法不能解决此类问题，使用的是 active-set 算法。如需更多帮助，请参阅文档中的选择算法。
@nasim，我用一些随机生成的数据检查了这段代码，所以我自己无法检查这个警告。正如我对这个问题的评论，最好发布你的矩阵。
我很乐意，但 A 矩阵不适合 cmets。我可以在哪里发帖？
@nasim，点击您问题下方的“编辑”链接，然后将数据复制到问题中。
@nasim，感谢您提供的固定数据。但是有了这些值，我得到了一个没有警告的正确解决方案。即使没有指定有限界，该解也是全正的；最小的 x 是 3.9e-5。

【解决方案3】：

当您有 x>0 条件时，这将是一个线性优化问题。解决这个问题的最佳算法是单纯形算法。这个想法是每个线性方程 aix=bi 提供一条线，这些线的组合提供一个多边形/多面体，答案是这个多面体/多边形的顶点之一。单纯形算法非常标准，有许多可用的函数和库可以计算。

【讨论】：