Create this particular SPARSE matrix to solve a quadratic program

Hello everyone,

I am trying to solve the following large scale quadratic program:

min (1/2*x'*Q*x + c'*x) subject to x >= 0

where the vector x is of dimension (n*d^2) where n is very large and d is small (in my problem: n = 220512 and d = 16).

The matrix Q is sparse, and is constructed as follows.

    one = ones(d,1);
    A = kron(eye(d),one');
    B = repmat(diag(one),1,d);
    M = A'*A + B'*B;

and Q is a *block diagonal* matrix where there are n identical blocks (in the diagonal), each block is the matrix M .

Now I would like to create Q as a sparse matrix (otherwise there will be a MEMORY problem when solving the original quadratic program, using quadprog for example).

Someone suggested me the following solution:

S = d^2*n;
Q = sparse(S, S); %'// preallocates with zeros
for b = 0:d^2:S-1
    Q(b+(1:d^2),b+(1:d^2)) = M;
end

However I think this solution is only suitable for small scale. I executed the code for n = 220512 and d = 16 an hour ago but it still keeps running (while the specs of my computer are not so bad: 16GB RAM, Intel Core i7-4770 3.4 GHz, quad-core).

Thank you in advance for any suggestions.

Answer 1

[Titus Edelhofer]

Hi,

I guess this will not work. If I take a look at M for d=16, it's 256*256 and requires (converted to sparse) about 0.12MByte. Just storing this 0.12MByte 220512 times would require 26 GByte of memory.

Nevertheless, here is what you would have to do:

MSparse = sparse(M);
% replicate: please don't test with n=220512 directly
Mblk = repmat({MSparse}, 1, n);
% convert do block diagonal
Q = blkdiag(Mblk{:});

Titus

[f10w]

Thanks Titus. Your solution works for me. The matrix is created and stored after about 15 minutes.

However, the matrix takes ~9GB of memory and thus slows down everything (although it takes only 1.1Gb when saved to file), so I guess there is no much hope to solve the original quadratic problem using quadprog (an example is given here).