Allocate a big matrix

I want to allocate a (10^5,10^4) matrix and after I want to fill the matrix in a for-loop with values. My problem is, that zeros(10^5,10^4) gives the error "matrix exceeds maximum array size preference." and sparse(10^5,10^4) makes the for-loop to slow. Has someone of you a solution of this problem?

Thanks a lot

[KSSV]

What is the purpose of allocating a matrix?

[Louis Wyss]

If I do not allocate the matrix, the for loop is also very slow. Matlab says I should allocate the matrix before the for-loop. I use the matrix for image reconstruction

[James Tursa]

You should not be filling in a sparse matrix in a for-loop. The usual technique advised is to save the values and indexes in separate variables and then combine into a sparse matrix at the end. Can you show us what your for-loop looks like?

[Louis Wyss]

numTrans = 128;
numWinkel = 720;
numPixX = 100;
numPixZ = 100;
  yT2 = zeros(numTrans*numWinkel,numPixZ*numPixX);
for k = 1:length(phi)
  for l = 1:numTrans
      for j = 1:numPixX
          for i = 1:numPixZ
              Int = A*exp(-(abs(xTransducer(l) - xiTransKoor(i,j,k))^2)/(a^2));
              if Int < 10^(-50)
                  yT2((k-1)*numTrans + l,(j-1)*numPixX + i) = 0;
              else
                  yT2((k-1)*numTrans + l,(j-1)*numPixX + i) = Int; 
              end
          end
      end
  end
end

Answer 1

[Walter Roberson]

Use spalloc() specifying the number of expected non-zero values.

Sparse matrices are more efficiently filled working down columns than across rows.

10^5 by 10^4 is 10^9, just under 1 giga-entries. If you were to use the default double() datatype, that would require 8 gigabytes to store. Your preferences have been configured to not allow you to store matrices that large. If you have enough memory for this (and for the other arrays you are using) then you could change your preferences. And if you can get away with using something other than double precision, then create the array that way. For example,

A = zeros(10^5, 10^4, 'int16');

would take only 1/4 of the storage.

[Louis Wyss]

Can I initialise the matrix as int16 and after fill up with double values

[Walter Roberson]

If you initialize as int16 and write double values then the double values will be converted to integers.

If you need to store double, then you cannot use this technique.

Are you trying to do something like subpixel resolution? If so then you simply do not have enough memory to work at that resolution.

[Louis Wyss]

I think it is possible a memory problem. Maybe I should try it on a PC with more memory.

Answer 2

[Jan]

You can increase the accepted size of allocated arrays in Matlab's preferences. If you have enough RAM to allocate an array of 8GB, there is no reason to restrict the arrays to smaller sizes.

Did you try to allocate the sparse matrix with a matching number of non-zero elements? Are you sure that the low speed of the loop is a problem of the sparse array? Perhaps there are other reasons. Perhaps somebody has a suggestion for improvements, when you post the relevant part of the code.

[Louis Wyss]

numTrans = 128;
numWinkel = 720;
numPixX = 100;
numPixZ = 100;
  yT2 = zeros(numTrans*numWinkel,numPixZ*numPixX);
for k = 1:length(phi)
  for l = 1:numTrans
      for j = 1:numPixX
          for i = 1:numPixZ
              Int = A*exp(-(abs(xTransducer(l) - xiTransKoor(i,j,k))^2)/(a^2));
              if Int < 10^(-50)
                  yT2((k-1)*numTrans + l,(j-1)*numPixX + i) = 0;
              else
                  yT2((k-1)*numTrans + l,(j-1)*numPixX + i) = Int; 
              end
          end
      end
  end
end
end

[Louis Wyss]

I tried it with sparse with smaller values for the loop and it was slower, but maybe you are right and my solution is very inefficient. I have not 8Gb Ram on my Laptop, so I have to search another solution. Maybe someone has a more efficient idea for my code. For loops are not the fastet but I was not able find a better solution.

[Walter Roberson]

What are some typical value for xTransducer(l) and xiTransKoor(i,j,k) ? So that we can test our ideas?

[Louis Wyss]

Typical values are single values. Most of the values are like 20.2355 or also integer. But I need the single format

[Walter Roberson]

Are any of your transducer values or xiTransKoor values complex valued? Because if not then the abs() part is not needed because of the ^2.

[Walter Roberson]

The vectorized form of your calculation appears to be:

yT2 = A*exp(-abs(bsxfun(@minus, xTransducer(:), reshape(xiTransKoor, [1 size(xiTransKoor)])).^2/a^2);
yT2 = reshape(yT2, numTrans*numWinkel, numPixZ*numPixX);

yT2(yT2 < 1E-50) = 0;

However that last step ends up building a temporary logical array the same size as the original, which could take a fair bit of memory. If you don't need the small values set to 0, then omit that last step. If you do need the small values set to 0, then you can reduce the temporary memory by using a for loop, something like

for K = 1 : size(yT2,2)
  yT2( yT2(:,K) < 1E-50, K) = 0;
end

[Jan]

Comments to your code:

• 10^(-50) is an expensive power operation while 1e-50 is a constant.
• Omit the line "yT2((k-1)*numTrans + l,(j-1)*numPixX + i) = 0;", because the elements are set to zeros already.
• Avoid calculating a^2 in each iteration by creating it as a constant before the loops.
• In

 Int = A*exp(-(abs(xTransducer(l) - xiTransKoor(i,j,k))^2)/(a^2));

you can omit the abs(), when the values of xTransducer and xiTransKoor are real.

• If you want to work with an 8GB array, prefer installing 16GB or more RAM, as a rule of thumb.
• I assume the bottleneck of your code is the expensive EXP() function. But you do not need this in many cases:

A*exp(-(abs(xTransducer(l) - xiTransKoor(i,j,k))^2)/(a^2)) < 1e-50
==> log(A) - (xTransducer(l) - xiTransKoor(i,j,k))^2/a^2 < log(1e-50)
==> abs(xTransducer(l) - xiTransKoor(i,j,k)) < sqrt((-log(1e-50) + log(A)) * a^2)

Then your code can be:

numTrans  = 128;
numWinkel = 720;
numPixX   = 100;
numPixZ   = 100;
yT2       = zeros(numTrans*numWinkel,numPixZ*numPixX);  % sparse!
a2        = a^2;
c         = sqrt((-log(1e-50) + log(A)) * a2);
for k = 1:length(phi)
  for l = 1:numTrans
    for j = 1:numPixX
      for i = 1:numPixZ
        d = abs(xTransducer(l) - xiTransKoor(i,j,k));
        if d >= c
          Int = A * exp(-d^2 / a2);
          yT2((k-1)*numTrans + l,(j-1)*numPixX + i) = Int;
        end
      end
    end
  end
end

Please check the calculations carefully, it is late in the night here. But if the output array conatins many zeros, reducing the number of expensive EXP() calls will accelerate the code remarkably.

[Louis Wyss]

Thanks for all your fast answers! I look them up and response if there are still some problems. I'm very thankful for all your help!!

[Louis Wyss]

Your recommedations made my program very faster, thanks again for your help. May I ask a second question about invert a matrix? I need to invert the yT2 matrix and I use the tikonov regularization. pinv() is not possible because the reconstruction is very bad when I use pinv(). I tried like in yT2 to split up the calculations but I want to ask you if you know a more efficient way

lambda = 250;  
tikhonovMatrix = lambda*speye(numTrans*numWinkel, numPixZ*numPixX);
tikhonovMatrixTrans = tikhonovMatrix';
tikhonovMatrixProd = tikhonovMatrixTrans*tikhonovMatrix;
yT2Trans = yT2';
yT2Prod = yT2Trans*yT2;
pseudoInverse_yT2 = yT2Prod + tikhonovMatrixProd;
bildVektor = reshape(newsinogram,[],1);
recBildVektor = yT2Trans*bildVektor;
recBildVektor2 = pseudoInverse_yT2\recBildVektor;

[Jan]

Please open a new thread for a new question.