Matrix version of quadprog?

9 views (last 30 days)
Eric Zhang
Eric Zhang on 30 Jun 2016
Edited: Eric Zhang on 1 Jul 2016
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?

Accepted Answer

Teja Muppirala
Teja Muppirala on 1 Jul 2016
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 Comment
Eric Zhang
Eric Zhang on 1 Jul 2016
Edited: Eric Zhang on 1 Jul 2016
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!

Sign in to comment.

More Answers (1)

John D'Errico
John D'Errico on 30 Jun 2016
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 Comment
Eric Zhang
Eric Zhang on 1 Jul 2016
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!

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!