How to speed up making large matrices, and reduce memory usage

How can I speed up making a large matrix, and also cut down memory usage?

I need to make a matrix which turns out to be in the order of 2,000 by 30,000 ish.

I can run it, but it'd be great to be able to speed it up, and also reduce memory use.

At the moment my code is as below, there length of C is 1800.

Creating Aineq is my final aim!

Aineq=[ones(length(C))]; %size of matrix to represent size one of the decision variables Q or G
Aineq1=[tril(Aineq),-tril(Aineq),tril(Aineq),-tril(Aineq),Aineq*0,Aineq*0,Aineq*0,zeros((length(C)),length(F))]; %1 matrix for first DV, then the next DV
Aineq2=[-tril(Aineq),tril(Aineq),-tril(Aineq),tril(Aineq),Aineq*0,Aineq*0,Aineq*0,zeros((length(C)),length(F))]; %2
Aineq3=[eye(length(C)),eye(length(C)),Aineq*0,Aineq*0,Aineq*0,Aineq*0,Aineq*0,zeros((length(C)),length(F))]; %3
Aineq4=[tril(Aineq,-1)*-1,tril(Aineq),tril(Aineq,-1)*-1,tril(Aineq,-1),Aineq*0,Aineq*0,Aineq*0,zeros((length(C)),length(F))]; %4
Aineq5=[Aineq*0,Aineq*0,Aineq*0,eye(length(C)),Aineq*0,Aineq*0,Aineq*0,zeros((length(C)),length(F))]; %5
Aineq6=[Aineq*0,Aineq*0,tril(Aineq,-1)*-1,tril(Aineq,-1),tril(Aineq,-1),tril(Aineq,-1)*-1,tril(Aineq,-1)*-1,zeros((length(C)),length(F))]; %6
Aineq7=[Aineq*0,Aineq*0,tril(Aineq)*-1,tril(Aineq),tril(Aineq),tril(Aineq)*-1,tril(Aineq)*-1,zeros((length(C)),length(F))]; %7
Aineq8=[Aineq*0,Aineq*0,tril(Aineq,-1),tril(Aineq,-1)*-1,tril(Aineq,-1)*-1,tril(Aineq,-1),tril(Aineq,-1),zeros((length(C)),length(F))]; %8
Aineq9=[Aineq*0,Aineq*0,tril(Aineq),tril(Aineq)*-1,tril(Aineq)*-1,tril(Aineq),tril(Aineq),zeros((length(C)),length(F))];%9
Aineq10=[Aineq*0,Aineq*0,eye(length(C)),Aineq*0,Aineq*0,Aineq*0,Aineq*0,zeros((length(C)),length(F))]; %10
Aineq11=[Aineq*0,Aineq*0,Aineq*0,eye(length(C)),Aineq*0,Aineq*0,Aineq*0,zeros((length(C)),length(F))]; %11
Aineq12=[Aineq*0,Aineq*0,Aineq*0,Aineq*0,Aineq*0,eye(length(C)),Aineq*0,zeros((length(C)),length(F))]; %12
Aineq13=[Aineq*0,Aineq*0,Aineq*0,Aineq*0,eye(length(C)),Aineq*0,Aineq*0,zeros((length(C)),length(F))]; %13
Aineq14=[Aineq*0,Aineq*0,Aineq*0,Aineq*0,Aineq*0,Aineq*0,eye(length(C)),zeros((length(C)),length(F))]; %14  
Aineq15=[Aineq*0,Aineq*0,Aineq*0,Aineq*0,Aineq*0,Aineq*0,eye(length(C))*-1,zeros((length(C)),length(F))]; %15


clear Aineq
% Aineq=[Aineq1;Aineq2;Aineq3;Aineq4;Aineq5;Aineq6;Aineq7;Aineq8;Aineq9;Aineq10;Aineq11;Aineq12;Aineq13;Aineq14;Aineq15];
Aineq=[Aineq1]; clear Aineq1
Aineq=[Aineq;Aineq2]; clear Aineq2
Aineq=[Aineq;Aineq3]; clear Aineq3
Aineq=[Aineq;Aineq4]; clear Aineq4
Aineq=[Aineq;Aineq5]; clear Aineq5
Aineq=[Aineq;Aineq6]; clear Aineq6
Aineq=[Aineq;Aineq7]; clear Aineq7
Aineq=[Aineq;Aineq8]; clear Aineq8 
Aineq=[Aineq;Aineq9]; clear Aineq9
Aineq=[Aineq;Aineq10]; clear Aineq10
Aineq=[Aineq;Aineq11]; clear Aineq11
Aineq=[Aineq;Aineq12]; clear Aineq12
Aineq=[Aineq;Aineq13]; clear Aineq13
Aineq=[Aineq;Aineq14]; clear Aineq14
Aineq=[Aineq;Aineq15]; clear Aineq15

Answer 1

[Bjorn Gustavsson]

This looks like a very sparse matrix, right? Then I would first figure out how to assign the blocks of Aineq from its components, something like this:

blksz = [3,13]; % or whatever they become
nBlks = 15;
Aineq = zeros(blksz*nBlks);
for i1 = 1:nBlks
  for i2 = 1:nBlks
    % Some logics here for selecting which block to fill Aineq with
    Aineq((1:blksz(1))+blksz*(i1-1),(1:blksz(2))+blksz(2)*(i2-1)) = CurrBlk;
  end
end

Once you've figured that one out, you should look at what components of every block is non-zero and start to build sparse blocks, and then concatenate them instead, that would be much more memmory efficient. Then the last step (or preferably the firts, but I find it easier to move to sparse representations of comlex matrices step-by-step), is to calculate vectors with the row and column-indices and the corresponding matrix element and then simply concatenate those into three arrays i_rows, i_cols, val and then call sparse:

Asparse = sparse(i_rows,i_cols,vals,nRows,nCols);

HTH

[In-chan Kim]

Thanks for your response!

So this would be in some ways the reverse of what I'm doing right?

Instead of making each segment of the overall matrix and putting them together, making the full matrix of zeros, then filling in the bits?

I'm unfamiliar with building sparse blocks and concatenating them.

I've had a search, and it seems like it "sequeezes out" the zeros. Just to understand correctly, does that mean literally getting rid of the zeros? Or are the zeros there, just in a different format which allows saving storage?

Because I would need those zeros there to do matrix multiplication.

Thanks!

[Bjorn Gustavsson]

1, Yeah, for step/stage 1, yes. You build all the submatrices at once and then start to concatenate them, at each concatenation-step matlab will have to re-allocate memmory - which is costly. It should be better to generate each block when it is needed and assign that block to where you want it in your Aineq matrix.

2, Yeah you can make a sparse matrix by squeezing out the zeroes. However for you you should go the other way around:

i_row = [1 1 2 2 3 4 6];
i_col = [1 3 2 4 1 5 2];
vals =  [3 1 4 1 5 9 2];
M1 = sparse(i_row,i_col,vals,6,6);
i_r = [1 3 5 2 4 6];
i_c = [6 5 2 1 4 3];
v =   [2 7 1 8 2 8];
M2 = sparse(i_r,i_c,v);
M = [M1;M2]
P = M1*M2
Mc = sparse([i_row,i_r+6],[i_col,i_c],[vals,v]);

So you generate the index and value arrays and then create the sparse matrix.

You really don't "need the zeroes" - matlab has a very efficient set of functions for operating on sparse arrays, so that will automatically be handled. To me it still takes too much thinking (for my comfort) to generate these arrays with ease, but when the matrices becomes big this is way way way more efficient than using full matrices.

If you need to solve equations with sparse matrices you should also have a look at the Factorize toolbox by Tim Davis, it is very good.

HTH

[Stephen23]

"Or are the zeros there, just in a different format which allows saving storage?"

The zeros are not stored explicitly in sparse classes, which saves on memory (for appropriate data).

"Because I would need those zeros there to do matrix multiplication."

No, you would not.

You can multiply the sparse matrix and MATLAB will return the correct result, that is the whole point of them (a sparse matrix that could not have basic numeric operations applied to it would be completely useless to anyone).

Create a small example and try it.

As Bjorn Gustavsson wrote, creating a large sparse matrix by defining all elements and then removing the zeros would be very inefficient. The efficient approach is to just define the non-zero elements and their indices, and call sparse once with those.

[Steven Lord]

To illustrate the savings you can get by storing an appropriate matrix as a sparse matrix, let's take a pretty favorable example for sparse: a matrix with only a few non-zero elements, the identity. We can build a full identity with the eye function and a sparse identity with speye.

>> A = eye(20);
>> B = speye(20);

The two matrices are equal

>> isequal(A, B)
ans =
  logical
   1

How much memory does MATLAB take to store them?

>> whos A B
  Name       Size            Bytes  Class     Attributes

  A         20x20             3200  double              
  B         20x20              488  double    sparse    

On the other hand, matrices with many non-zero elements can actually take more memory to store as sparse. The worse case scenario is a sparse with all non-zero elements.

>> C = magic(20);
>> D = sparse(C);
>> whos C D
  Name       Size            Bytes  Class     Attributes

  C         20x20             3200  double              
  D         20x20             6568  double    sparse    

[Bjorn Gustavsson]

...(to continue the efficiency illustrations) the speedup of matrix-multiplications are not to be fogotten:

R = sprandn(3000,3000,0.005);
fR = full(R);
tic,M2 = fR*fR;toc
% gave me: Elapsed time is 0.397631 seconds.
tic,M2 = R*R;toc  
% and: Elapsed time is 0.048307 seconds.

The speed-up-factor depends on the sparsity-factor of the matrix but can become substantial if the matrix is "very" sparse.

[In-chan Kim]

Brilliant! Thanks for everyone's responses! Trying this now!

[In-chan Kim]

Thanks everyone, real good outputs! About 3 times less memory and faster!

Name          Size                    Bytes  Class     Attributes

  A         27003x12603            2722550472  double              
  spr       27003x12603             829912320  double    sparse    
Elapsed time is 4.461101 seconds.
Elapsed time is 15.440768 seconds.

[In-chan Kim]

Is there any way to see how much memory is being used by a code being run? And which parts are taking up a lot of memory?

I know with profiler you can do that for the time it takes.

Also, I thought I'd share what code I went with

Aineqinit=zeros(15*length(C)+length(F),7*length(C)+length(F)); % full size of Aineq based on the make up of the 16 constraints
Aineqbase=[ones(length(C))];%size of matrix to represent size one of the decision variables Q or G
Aineqinit(1:length(C),1:7*length(C)+length(F))=[tril(Aineqbase),-tril(Aineqbase),tril(Aineqbase),-tril(Aineqbase),Aineqbase*0,Aineqbase*0,Aineqbase*0,zeros((length(C)),length(F))]; %1 matrix for first DV, then the next DV
Aineqinit(1*length(C)+1:2*length(C),1:7*length(C)+length(F))=[-tril(Aineqbase),tril(Aineqbase),-tril(Aineqbase),tril(Aineqbase),Aineqbase*0,Aineqbase*0,Aineqbase*0,zeros((length(C)),length(F))]; %2
Aineqinit(2*length(C)+1:3*length(C),1:7*length(C)+length(F))=[eye(length(C)),eye(length(C)),Aineqbase*0,Aineqbase*0,Aineqbase*0,Aineqbase*0,Aineqbase*0,zeros((length(C)),length(F))]; %3
Aineqinit(3*length(C)+1:4*length(C),1:7*length(C)+length(F))=[tril(Aineqbase,-1)*-1,tril(Aineqbase),tril(Aineqbase,-1)*-1,tril(Aineqbase,-1),Aineqbase*0,Aineqbase*0,Aineqbase*0,zeros((length(C)),length(F))]; %4
Aineqinit(4*length(C)+1:5*length(C),1:7*length(C)+length(F))=[Aineqbase*0,Aineqbase*0,Aineqbase*0,eye(length(C)),Aineqbase*0,Aineqbase*0,Aineqbase*0,zeros((length(C)),length(F))]; %5
Aineqinit(5*length(C)+1:6*length(C),1:7*length(C)+length(F))=[Aineqbase*0,Aineqbase*0,tril(Aineqbase,-1)*-1,tril(Aineqbase,-1),tril(Aineqbase,-1),tril(Aineqbase,-1)*-1,tril(Aineqbase,-1)*-1,zeros((length(C)),length(F))]; %6
Aineqinit(6*length(C)+1:7*length(C),1:7*length(C)+length(F))=[Aineqbase*0,Aineqbase*0,tril(Aineqbase)*-1,tril(Aineqbase),tril(Aineqbase),tril(Aineqbase)*-1,tril(Aineqbase)*-1,zeros((length(C)),length(F))]; %7
Aineqinit(7*length(C)+1:8*length(C),1:7*length(C)+length(F))=[Aineqbase*0,Aineqbase*0,tril(Aineqbase,-1),tril(Aineqbase,-1)*-1,tril(Aineqbase,-1)*-1,tril(Aineqbase,-1),tril(Aineqbase,-1),zeros((length(C)),length(F))]; %8
Aineqinit(8*length(C)+1:9*length(C),1:7*length(C)+length(F))=[Aineqbase*0,Aineqbase*0,tril(Aineqbase),tril(Aineqbase)*-1,tril(Aineqbase)*-1,tril(Aineqbase),tril(Aineqbase),zeros((length(C)),length(F))];%9
Aineqinit(9*length(C)+1:10*length(C),1:7*length(C)+length(F))=[Aineqbase*0,Aineqbase*0,eye(length(C)),Aineqbase*0,Aineqbase*0,Aineqbase*0,Aineqbase*0,zeros((length(C)),length(F))]; %10
Aineqinit(10*length(C)+1:11*length(C),1:7*length(C)+length(F))=[Aineqbase*0,Aineqbase*0,Aineqbase*0,eye(length(C)),Aineqbase*0,Aineqbase*0,Aineqbase*0,zeros((length(C)),length(F))]; %11
Aineqinit(11*length(C)+1:12*length(C),1:7*length(C)+length(F))=[Aineqbase*0,Aineqbase*0,Aineqbase*0,Aineqbase*0,Aineqbase*0,eye(length(C)),Aineqbase*0,zeros((length(C)),length(F))]; %12
Aineqinit(12*length(C)+1:13*length(C),1:7*length(C)+length(F))=[Aineqbase*0,Aineqbase*0,Aineqbase*0,Aineqbase*0,eye(length(C)),Aineqbase*0,Aineqbase*0,zeros((length(C)),length(F))]; %13
Aineqinit(13*length(C)+1:14*length(C),1:7*length(C)+length(F))=[Aineqbase*0,Aineqbase*0,Aineqbase*0,Aineqbase*0,Aineqbase*0,Aineqbase*0,eye(length(C)),zeros((length(C)),length(F))]; %14  
Aineqinit(14*length(C)+1:15*length(C),1:7*length(C)+length(F))=[Aineqbase*0,Aineqbase*0,Aineqbase*0,Aineqbase*0,Aineqbase*0,Aineqbase*0,eye(length(C))*-1,zeros((length(C)),length(F))]; %15
% formulate LHS of constraint 16
u = zeros((numel(stcounter)+1)/2, sum(stcounter));
for i = 1:2:numel(stcounter)
    u(((i+1)/2), :) = repelem([0, 1, 0], [sum(stcounter(1:i-1)), stcounter(i), sum(stcounter(i+1:end))]);
end
Aineqinit(15*length(C)+1:15*length(C)+length(F),1:7*length(C)+length(F))=[u,zeros(length(F),length(C)),zeros(length(F),length(C)),zeros(length(F),length(C)),zeros(length(F),length(C)),zeros(length(F),length(C)),zeros(length(F),length(C)),eye(length(F))*-1*dtih.*PEVmaxstay.*staycounter]; %16
Aineq=sparse(Aineqinit);