Matrix version of quadprog?

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?

Answer 1

[Teja Muppirala]

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.

[Eric Zhang]

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!

Answer 2

[John D'Errico]

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.

[Eric Zhang]

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!