Matrix version of quadprog?
4 views (last 30 days)
Show older comments
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?
0 Comments
Accepted Answer
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.
More Answers (1)
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.
See Also
Categories
Find more on Linear Least Squares 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!