Matrix version of quadprog?

10 ビュー (過去 30 日間)
Eric Zhang
Eric Zhang 2016 年 6 月 30 日
編集済み: Eric Zhang 2016 年 7 月 1 日
I have a standard quadratic programming problem with equality constraint as outlined here, except that instead of vector x, I am optimizing over matrix X, whose each row is a data point. The objective function is tr(X^THX).
quadprog seems to be unable to handle matrix X. Am I missing something, or I have to go find some unofficial packages? If so, can one point me to such a solver?

採用された回答

Teja Muppirala
Teja Muppirala 2016 年 7 月 1 日
You can convert this into a form usable by QUADPROG by noting the identity:
tr(X'*H*X) == vec(X)' * kron(I, H) * vec(X)
where I'm defining vec(X) = X(:) (basically reshape the n x n matrix into a vector of n^2 elements.)
An example will make this clear.
----------------------------------------------
Problem:
Minimize tr(X'*H*X)
subject to
Aeq*x = beq
where x is defined as X(:), or equivalently X = reshape(x,size(H))
----------------------------------------------
Solution:
%%Step 0. First make some random data
rng(0);
% Symmetric positive semidefinite matrix H
n = 5; % H is size n x n
H = randn(n,n);
H = H'*H;
% Make some constraints
nC = 3; % Number of constraints
Aeq = randn(nC,numel(H));
beq = randn(nC,1);
%%Method 1: General constrained optimization solver FMINCON
unvec = @(x) reshape(x,size(H));
traceFun = @(x) trace( unvec(x)'*H*unvec(x) );
[x_opt,f_opt] = fmincon(traceFun, zeros(numel(H),1),[],[],Aeq,beq);
%%Method 2: QUADPROG
% These two lines are equivalent to Q = kron(eye(size(H)),H);
% Use sparse to save space.
Q = repmat({H},size(H,1),1);
Q = spblkdiag(Q{:});
% Use 2*Q because QUADPROG multiplies by 1/2
[x_optQ,f_optQ] = quadprog(2*Q,[],[],[],Aeq,beq);
%%3. Verify that they both give the same answers
f_opt
f_optQ
unvec(x_opt)
unvec(x_optQ)
When you run, this you see that you get the same answer either way:
f_opt =
0.0527
f_optQ =
0.0527
ans =
-0.0832 0.0434 0.0059 -0.0198 -0.0223
0.0795 -0.0289 -0.0218 0.0385 0.0260
0.0948 -0.0718 0.0457 -0.0320 0.0092
-0.2815 0.1671 -0.0514 -0.0101 -0.0543
0.1309 -0.1191 0.1142 -0.0878 0.0009
ans =
-0.0832 0.0434 0.0059 -0.0198 -0.0223
0.0795 -0.0289 -0.0218 0.0385 0.0260
0.0948 -0.0718 0.0457 -0.0320 0.0092
-0.2815 0.1671 -0.0514 -0.0101 -0.0543
0.1309 -0.1191 0.1142 -0.0878 0.0009
QUADPROG will converge much faster than FMINCON, and should work for moderate problem sizes (n ~ 100), but for much larger problem sizes, making the kronecker product (even when sparse) will result in an out of memory error.
  1 件のコメント
Eric Zhang
Eric Zhang 2016 年 7 月 1 日
編集済み: Eric Zhang 2016 年 7 月 1 日
This is eye-opening. Hats off! This teaches me to think more about massaging my objective function into one suitable for existing MATLAB functions. Thanks a lot!

サインインしてコメントする。

その他の回答 (1 件)

John D'Errico
John D'Errico 2016 年 6 月 30 日
This problem is not a form that quadprog can solve. You can write it in a more computationally efficient way though as an objective function, perhaps something that fmincon could use.
f = sum(sum(X.*(H*X),1));
This computes ONLY the diagonal elements that are needed in the trace, then sums them. For example:
X = randn(1000,1000);
H = rand(1000,1000);
timeit(@() trace(X'*H*X))
ans =
0.10347
timeit(@() sum(sum(X.*(H*X),1)))
ans =
0.053614
trace(X'*H*X)
ans =
5.0902e+05
sum(sum(X.*(H*X),1))
ans =
5.0902e+05
As you can see, the two results are the same, but one is a bit faster than the other.
  1 件のコメント
Eric Zhang
Eric Zhang 2016 年 7 月 1 日
Thanks John for the answer! Teja's answer managed to reduce my problem to one that quadprog can solve, so I accepted his answer. However, your answer indeed shed light on my another problem at hand. So thanks a lot!

サインインしてコメントする。

カテゴリ

Find more on Linear Programming and Mixed-Integer Linear Programming in Help Center and File Exchange

タグ

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by