使用 fmincon 在 matlab 上优化矩阵

Question

我有一个 30x30 矩阵作为基础矩阵 (OD_b1)，我还有两个基础向量（bg 和 Ag）。我的目标是优化尺寸为 30X30 的矩阵 (X)，使得： 1）向量（bg）与所有列的和向量之间的平方差被最小化。 2）向量（Ag）与所有行之和向量之间的平方差被最小化。 3)矩阵(X)和矩阵(OD_b1)的元素的平方差被最小化。

方程的数学形式如下：

我试过这个：

 fun=@(X)transpose(bg-sum(X,2))*(bg-sum(X,2))+ (Ag-sum(X,1))*transpose(Ag-sum(X,1))+sumsqr(X_b-X);
[val,X]=fmincon(fun,OD_b1,AA,BB,Aeq,beq,LB,UB)

我没有收到错误，但它似乎卡住了。

是我变量太多还是有别的原因？

提前致谢

Answer 1

这是一个简单的、无约束的最小二乘问题，因此有一个简单的解，可以表示为线性系统的解。

我将向您展示 (1) 解决此问题的精确有效的方法以及 (2) 如何使用 fmincon 解决。

精确、高效的解决方案：

问题设置

为了让我们在同一个页面上，我将变量初始化如下：

n = 30;
Ag  = randn(n, 1);    % observe the dimensions
X_b = randn(n, n);
bg  = randn(n, 1);

代码：

A1 = kron(ones(1,n), eye(n));
A2 = kron(eye(n), ones(1,n));
A = (A1'*A1 + A2'*A2 + eye(n^2));
b = A1'*bg + A2'*Ag + X_b(:);
x = A \ b;                      % solves A*x = b
Xstar = reshape(x, n, n);

为什么有效：

我首先重新表述了您的问题，因此目标是向量x，而不是矩阵X。观察z = bg - sum(X,2) 等价于：

x = X(:)                        % vectorize X
A1 = kron(ones(1,n), eye(n));   % creates a special matrix that sums up
                                % stuff appropriately
z = A1*x;

同样，A2 的设置使A2*x 等同于Ag'-sum(X,1)。那么你的问题就相当于：

minimize (over x) (bg - A1*x)'*(bg - A1*x) + (Ag - A2*x)'*(Ag - A2*x) + (y - x)'*(y-x) 其中y = Xb(:)。也就是说，y 是 Xb 的矢量化版本。

这个问题是凸的，一阶条件是最优的充要条件。取关于 x 的导数，该方程将定义您的解决方案！几乎等效的示例数学示例（但稍微简单的问题如下）：

minimize(over x) (b - A*x)'*(b - A*x) + (y - x)' * (y - x)

重写目标：

b'b- b'Ax - x'A'b + x'A'Ax +y'y - 2y'x+x'x

相当于：

minimize(over x) (-2 b'A - 2y'*I) x + x' ( A'A + I) * x

一阶条件为：

(A'A+I+(A'A+I)')x -2A'b-2I'y = 0
(A'A+I) x = A'b+I'y

你的问题本质上是一样的。它具有一阶条件：

(A1'*A1 + A2'*A2 + I)*x = A1'*bg + A2'*Ag + y

如何用 fmincon 解决

您可以执行以下操作：

f = @(X) transpose(bg-sum(X,2))*(bg-sum(X,2)) + (Ag'-sum(X,1))*transpose(Ag'-sum(X,1))+sum(sum((X_b-X).^2));
o = optimoptions('fmincon');%MaxFunEvals',30000);
o.MaxFunEvals = 30000;
Xstar2 = fmincon(f,zeros(n,n),[],[],[],[],[],[],[],o);

然后您可以使用以下方法检查答案是否相同：

normdif = norm(Xstar - Xstar2)

您可以看到差距很小，但基于线性代数的解决方案更精确：

gap = f(Xstar2) - f(Xstar) 

如果 fmincon 方法挂起，请尝试使用较小的 n 以确信我的基于线性代数的解决方案更精确、更快等... n = 30 正在解决 30^2 = 900 变量优化问题：不容易。使用线性代数方法，您可以达到 n = 100（即 10000 个变量问题）甚至更大。

Answer 2

我可能会使用quadprog 使用以下重新表述将其作为 QP 解决（保持目标尽可能简单以使问题“减少非线性”）：

min sum(i,v(i)^2)+sum(i,w(i)^2)+sum((i,j),z(i,j)^2)
v = bg - sum(c,x)
w = ag - sum(r,x)
Z = xbase-x

QP 求解器更精确（没有使用有限差分的梯度）。这种方法还允许您添加额外的边界以及线性等式和不等式约束。

另一个明确形成一阶条件的建议也是一个很好的建议：它也没有不精确梯度的问题（一阶条件是线性的）。由于它的灵活性，我通常更喜欢二次模型。