How to fill up a huge/very large matrix, with elements from another matrix in a time efficient way?

Hi there,

I am doing some calculations for some graph/network with vertices that have certain connections/edges to each other.

My goal is to set-up a very large matrices with dimensions ranging from roughly 90000x300(works but already takes up ~4.6 seconds to create) up to 9,000,000 (9mln)x30000 elements in an efficient way, such that is does not take forever to fill up the matrix. Most elements are, however, zeros so I tried using sparse matrices, but sparse matrices are not suitable if you want to fill them up element-wise.

The size of the matrices depend on a discretization parameter(delta) that I choose, I need to make it as small as possible, down to 10^-4. My function is now as follows(I use Hungarian notation):

function mZ_j = createZ_j(iN, iP, iD_j, mX, vPA)
% Function that creates matrix Z_j. 
% Dimensions: (n-p)x(p*d_j)+1, d_j can go up to 6, but is 2 on average

% input/starting parametes are:
% delta = 0.0001, s=0.5, T=900(15 min in seconds), d_j=2

iRows = iN-iP; % n:= T/delta, p:=s/delta
iCol  = (iP*iD_j);
mZ_j = zeros(iRows,iCol); % (n-p)x(p*d_j) matrix
mX_j = mX(:,vPA);        % matrix were elements will be taken from, vPA(size iD_j) has the column indexes
p = iP+1;                % indicator for each row; it is an autoregressive process and each row you regress on p lags, but every row you take one next step in time
for i=1:iRows
    for j=1:iP
        iCStart = ((j-1)*iD_j)+1; iCEnd = j*iD_j;
        mZ_j(i,iCStart:iCEnd) = mX_j(p-j,:);
    end
    p = p+1;
end

Does anyone see or know how I can assign the elements from mX_j more efficiently to mZ_j? The operations themselve do not take a lot time, it is just the amount of iterations that it has to do.

Thanks!

[Walter Roberson]

https://www.mathworks.com/help/matlab/math/accessing-sparse-matrices.html#br7572b

[Guillaume]

Questions:

• Is mX sparse? Otherwise, I don't see where the sparsity comes from in mZ_j.
• isn''t iD_j equal to numel(vPA)
• You're filling up your result and reading your input by rows. It would be faster to do it by columns. could the mX input be generated transposed and the mz_J output used transposed? Even better would be to have mx = fliplr(mx.')

[Jeroen]

• mX is not defined as a sparse matrix but IS sparse. I am reluctant to use sparse matrices, because later on I need to do matrix multiplication and taking inverses. When the matrices are sparse, these operations go way slower--> maybe you know a solution to this? I need to do: later on.
• iD_j is the amount of 'parent' nodes/vertices, and yes this is equal to the size of vPA. vPa contains the indices of the columns in mX that correspond to these ' parent' nodes.
• the output of mZ_j transposed makes all the next operations outside this function invalid

[Guillaume]

I forgot to ask:

• is iN equal to size(mX, 1)+1 or in any way related to the height of mX

It's fairly trivial to construct mz_J without a loop as a full matrix with simple indexing. The most expensive operation, if you're required to operate along the rows, is going to be some transpose to temporary operate along columns. However, as point out by Steven, a 9 million by 30000 matrix is going to require massive amounts of memory.

As for the math question, I'm the wrong person to ask I'm afraid.

[Jeroen]

• iN~T/delta=9000000(9m) and mX is iNx6. So yes they are related

Yes, I understand that it will take up a lot of memory and that I therefore need to use sparse matrices. But I do not now how to, I never worked with sparse matrices

[Walter Roberson]

Inverse of a sparse matrix is typically dense. Inverse of a nonsquare matrix is not defined. But perhaps you are doing something like A*A' that will be square but possibly large.

Inverse of a large matrix that is not just 0 and 1 will typically have terrible condition number and give numeric garbage. Don't use inv unless you are forced to by external requirements. Use other techniques such as \ instead.

Answer 1

[Steven Lord]

Your 90000-by-300 matrix isn't too big. If it's a full real double precision matrix, it'll take up a couple hundred MB of contiguous space.

Your 9 million by 30000 matrix better be sparse or you will need a computer with a very large amount of memory. One copy of a full real matrix that size will require just shy of 2 TB of contiguous space.

Rather than filling the sparse matrix one element at a time, build three vectors: one of row indices, one of column indices, and one of data values. Call sparse once at the very end to assemble the matrix from those three vectors. See the "Sparse Matrix of Nonzeros with Specified Size" example on the documentation page for the sparse function for an illustration of the technique. The "Accumulate Values into Sparse Matrix" example may also be of interest.

[Jeroen]

Thank you.

I do not know in advance which indices get zero or nonzero values, and what these nonzero values are. This all depends on mX_j which changes on every run of the program. I just take 'blocks' of size 1xp*d_j from mX_j, fill it in mZ_j and go to the next row. The order does matter because it is an autoregressive process. Could you maybe explain how you would tackle this?

[Walter Roberson]

If there are a different number per block then preallocate cells by number of blocks and gather the information specific to the block. Later horzcat the parts together into vector that are then passed to sparse.

[Jeroen]

There is a pattern in the matrix. The first row minus the elements in the last d_j columns is shifted to the second row and shifted one place to the righte, e.g., row(1,end-1)-->row(2,2:end). Then only a new value is added in the first place of the new row. So if d_j=2 I need the matrix(forget about the vector of ones in the last column, it will be added at the very last):

I try to understand what you mean with how to solve it, but I cannot seem to find the way to do it..

[Jeroen]

I now have the following code, it works a lot faster than the previous one, but still it will not work for 9mlnx30000 matrices:

function mZ_j = createZ_j_new(iN, iP, iD_j, mX, vPA)
% Function that creates  matrix Z_j
% Dimensions: (n-p)x(p*d_j)+1 
% delta=0.0001, s=0.5,T=900(15 min in seconds), d_j is on average 2 but can go up to 6, vPa contains indices of columns in mX that are relevant.
% n:=T/delta=9000000, p:=s/delta=5000, d_j=6 (maximum)
% createx indexes to select from mX
newVec=(iP:1:iN-1)';
newRow = 0:1:(iP-1);
newRow = repmat(newRow,iD_j,1);
newRow = newRow(:)';
newMat = repmat(newVec,1,(iP*iD_j));
newMat2=bsxfun(@minus,newMat,newRow);
mIndicesJ=repmat(vPA',iN-iP,iP);
mZ_j=mX(sub2ind(size(mX),newMat2,mIndicesJ));
iRows = iN-iP;
mZ_j = [mZ_j ones(iRows,1)];    % add column of ones to Z_j
end

I am new to sparse matrices and I do not know how to translate this into sparse matrix operations. Can you help me?

[Steven Lord]

Here's a simple sparse example that you can use to try to determine how to build your sparse matrix. I'm going to build this matrix:

>> A = [1 0 0 0 6; 0 2 0 7 0; 0 0 11 0 0; 0 9 0 4 0; 10 0 0 0 5]

A =

     1     0     0     0     6
     0     2     0     7     0
     0     0    11     0     0
     0     9     0     4     0
    10     0     0     0     5

First, for the main diagonal the elements are in rows 1:5, columns 1:5. The values are 1:5. Yes, I know the center element is 11 not 3. Wait for the next step.

v = (1:5).';
mainDiagonal = [v, v, v];

For the antidiagonal, the elements are in rows 1:5, columns 5:-1:1. The values are 6:10. When you specify multiple values for the same row and column indices in your vector, sparse accumulates them. So the 3 in position (3, 3) from the main diagonal and the 8 in position (3, 3) from the antidiagonal will add up to 11. [It's almost like I designed the example that way ;)]

antiDiagonal = [v, flip(v), v+5];

Let's assemble those into a matrix whose first column is row indices, whose second is column indices, and whose third is values.

allElements = [mainDiagonal; antiDiagonal];

Finally, build and compare

S = sparse(allElements(:, 1), allElements(:, 2), allElements(:, 3))
full(S)

In this case, S is dense enough (9/25 or 36% of the elements are nonzero) that it actually takes up more memory than A does. If you were to scale this up to something a bit larger:

v2 = (1:99).';
S2 = sparse([v2; v2], [v2; flip(v2)], [v2; v2+99]);
A2 = full(S2);
whos A2 S2

You'll see you save a bunch of memory (and typing; specifying all the 0's in A2 like I did in A would have been quite tedious and prone to misplacement.) In constructing S2, I needed to build three full matrices instead of just the one that I would have needed to build A2. However, the full matrices I used to build S2 are each 198-by-1 (597 elements total) while the full matrix A2 is 99-by-99 (9801 elements.) I used an odd size so the center element accumulates values from both the main and anti diagonals to demonstrate that capability.

[Jeroen]

Thank you!