How to create a large matrix using another matrix

Hello everyone, I want to make a large matrix (10^7 x 10^7) but it needs to have the following matrix being repeated around its main diagonal.

a=4;
A=[1     0     a+1     0; 
  a+2    2      0     a+1; 
   0    a+2     3      0; 
   0     0     a+2     4];

As you can see it is a main diagonal with 1:10^7 and 1 row lower repeating the number 6 and 2 rows higher repeating the number 5.

Everything I have tried turns out to be huge in terms of memory and unable to be performed like that. I suppose the trick is by somehow making use of a sparse matrix, but I cannot get it to work properly.

Thanks a lot in advance!

[Steven Lord]

What are you trying to do with this large matrix? Are you trying to use it to solve a system of equations? If so consider instead using one of the iterative solvers in MATLAB. Most if not all of those solvers can operate on either an explicit coefficient matrix A or a function that performs one or both of A*x or A'*x, where x is a column vector. If your coefficient matrix has a pattern (as yours apparently does) you may be able to compute A*x and A'*x without explicitly forming A.

[Petros Tsitouras]

Thanks for your interest Steven!

The main goal is to solve for x in the following equation:

B*x=C

Where C=ones(10^7,1)

Finding the inv(B) to solve it is not possible at this size.

[Walter Roberson]

You should rarely be using inv() on any matrix that is more than about 4 x 4, and even then only on symbolic matrices. In other cases you should be using the \ operator. The \ operator should be able to detect that it is a band-restricted sparse system and use a more efficient solver.

[Petros Tsitouras]

Noted! But yet again x=C\B displays the same error concerning output's size.

[Walter Roberson]

What error is that?

a=4;
N=1E7;
Adiag=(1:N).';
A1=ones(N,1)*(a+2);
A2=ones(N,1)*(a+1);
Acom=[A2, Adiag, zeros(N,1), A1];
B=spdiags(Acom,-1:2,N,N);

C = ones(N,1);

sol = B\C;

On my system the solution is found about 1.6 seconds with no memory problems. Memory use peaks about less than 6 gigabytes on my system.

[Petros Tsitouras]

Found out the difference, I have calculated the Acom as follows:

Adiag=(1:N);
A1=ones(1,N)*(a+2);
A2=ones(1,N)*(a+1);
Acom=[A1; zeros(1,N); Adiag; A2];
B=spdiags(Acom',-1:2,N,N);

Which results to

Acom:

instead of:

and B:

instead of:

Not sure which approach is the mathematically correct one.

[Walter Roberson]

Your original request shows the a+2 below the diagonal, so anything that ends up with the 6 above the diagonal is a wrong approach ;-)

The approach I used of constructing columns instead of rows has the advantage of not needing to transpose Acom, and so is more efficient.

[Petros Tsitouras]

Then I will have to aggree with you and tell you a huge thanks once more!! Have a good day!

Answer 1

[Walter Roberson]

https://www.mathworks.com/help/matlab/ref/spdiags.html#mw_6de431d3-8887-4b16-b947-86a610ddea8c

[Petros Tsitouras]

Thanks for your answer Walter! It was really helpful

I tried the following:

a=4;
Adiag=1:10^3;
A1=ones(1,10^3)*(a+2);
A2=ones(1,10^3)*(a+1);
Acom=[A1; zeros(1,10^3); Adiag; A2];
B=spdiags(Acom',1,10^3,10^3);

The result is a 1000x1000 matrix but is way too big to display in matlab. So I tried at 100x100 and the result is this:

val =

   (1,2)        6
   (2,3)        6
   (3,4)        6
   (4,5)        6
   (5,6)        6
   (6,7)        6
   (7,8)        6
   (8,9)        6
   (9,10)       6
  (10,11)       6
  (11,12)       6
....... until (99,100)

[Walter Roberson]

your d should not be 1, it should be -1:2

[Petros Tsitouras]

Got it! Thanks a alot for the help! Is there anyway to upscale this to a 10^7 x 10^7 matrix?

[Walter Roberson]

a=4;
N=1E7;
Adiag=(1:N).';
A1=ones(N,1)*(a+2);
A2=ones(N,1)*(a+1);
Acom=[A2, Adiag, zeros(N,1), A1];
B=spdiags(Acom,-1:2,N,N);