# Block diagonal matrices in MATLAB?

• MATLAB

## Main Question or Discussion Point

So, using MATLAB, I'm trying to create the block diagonal matrix for the central difference approximation to the 2D laplacian operator in rectangular coordinates, and I've stumbled on to what looks like a pretty useful function. It's called "blkdiag", and it works by taking all of the input arguments, which are matrices, and placing them on the diagonal of a larger sparse matrix.

The problem I have with this function however, is that it doesn't create a block diagonal matrix of arbitrary size--it only concatenates the numer of arguments you input. In my case, all of the blocks are identical (viz., they are the tribanded matrix, B[1,-2,1]), and I want to place N-1 of these blocks on the diagonal of the larger matrix.

Obviously, I'm not going to manually type in the block N-1 times, so I was hoping that maybe someone out there would have a clever solution to this problem. If there are other ways of creating the matrix I'm looking for, that's great too, but it would be neat if I could get this particular method to work. I'm rather new to MATLAB so forgive me if this seems like a dumb question.

Related MATLAB, Maple, Mathematica, LaTeX News on Phys.org
Run a for loop and give it your B matrix. It will create another matrix (say C) with B repeated N-1 times and then give this C to blkdiag.

-- AI

If you have huge matrices, such a loop can be extremely long to run (order n^3, as you have to reallocate memory and copy at each step). As this page is the first to pop in google for "matlab block diagonal", this might be helpful to others:
Code:
>> a=cell(1,3)
a =
[]     []     []
>> [a{:}]=deal(sparse(rand(2)))
a =
[2x2 double]    [2x2 double]    [2x2 double]
>> b=blkdiag(a{:});
>> full(b)
ans =
0.3529    0.0099         0         0         0         0
0.8132    0.1389         0         0         0         0
0         0    0.3529    0.0099         0         0
0         0    0.8132    0.1389         0         0
0         0         0         0    0.3529    0.0099
0         0         0         0    0.8132    0.1389
The "sparse" is not required (but recommended for huge matrices). This also works if your blocks are not all the same, as long as you can have them in a cell array.

Pierre-André

Hi,
I have to write a block tridiagonal matrix of dimension 100 in matlab and I tried every thing but it is not working for me. I have a matrix B of dimension 10 which is tridiagonal and 4 on the main diagonal and super diagonal and -6 on the subdiagonal. From this B I have to develop a tridiagonal matrix A of dimension 100 whose leading diagonal is my above matrix B and whose upper and lower diagonals are Matrix -I. I will be grateful if anyone can help me how to write it in matlab and especially as a part of M-file.

Thanks
mfarooq52003

Suppose your matrix is B the block that you want to repeat 100 times, then

Code:
BL = kron(eye(100),B);
gives the result for NoOne.

mfarooq, try to use the following commands,

Code:
diag(-6*ones(1,14),-1) + diag(4*ones(1,15)) + diag(4*ones(1,14),1)
and modify according to your needs.

Hi trambolin

Thank you very much for your help and such a prompt response.

Regards

Farooq

Hi,

Now if I have a tridiagonal matrix B with main diagonal 4, upper diagonal 4 and lower diagonal -6 and I have to creat a matirx A with main diagonal the above matrix B, upper and lower diagonals Identity matrix I of the same dimension B, how we can creat that. My actual problem is writing this A in the above mentioned problem.

I will be grateful for your help.

Cheers
Farooq

Then you use the same trick over and over again.

to put the B matrices on the diagonal

Code:
BL = kron(eye(n),B);
then to put the identity matrices on the super on sub diagonals, let n is the number dictating how many times B is repeated

Code:
B = rand(3);
n = 10;
BL = kron(eye(n),B);
[p,q] = size(B);
BLI = BL + diag(ones(1,n*q-p),-p) + diag(1*ones(1,n*q-p),p);

Hi Trambolin,
Thank you once again for your help and quick response. It worked this time and my matrix A is now block tridiagonal. I just made the following changes according to my matrix structure,

BLI = BL - diag(ones(1,n*q-p),-p) - diag(1*ones(1,n*q-p),p); i.e changed the positive signd to negatives.

Thank you very much my friend once again.

Regards

Farooq

Hi,

i have blocktridiagonal matrix A of order n^2 * n^2.
all the entries are at the super, main and sub diagonal, that are tridiagonal matrices themselves, of order n*n.
i have stored all the entries(matrices) using 3D array.
Now i am lost, how to put them back in A as tridiagonals, as A is 2D matrix.
secondly, can i avoid 3D arrays. I have to pass A as argument in the triblocksolve algo by Thomas in Matlab.

Plz help, its urgent.

Thanks.
jave

Greetings ~

I am similarly new to Matlab & looking to define a matrix whose block-diagonal (2x2) matrices are identical, and whose off-diagonal elements are uniform. From the above I think I can sort out how to get the block-diagonal first, but I can't for the life of me sort how to set all the other elements without altering the (block)-diag. elements.

Any elegant solutions?!

Thank you.

Rax

First, my apologies to Jave I did not see his post and quite some time past.

Rax, what do you mean by uniform?

The matrix P represents probabilities associated with Ising spins.I think it is actually called a rank 4 tensor in 2D (?) - it has NxN 2x2 matrices as its elements that describe the probabilities of finding any given pair of spins aligned or anti-aligned.

As such, the diagonal 2x2 elements of P have a different probability distribution (viz. P(ii) = [1/2 0; 0 1/2] vs. P(ij) = [ 1/4 1/4 ; 1/4 1/4], for i neq. j). This is just for a uniform probability distribution, but ultimately I will be looking at more complicated distributions, so the more control over this entire matrix I have, the better.

At the moment, I've resolved the issue in the following way:

L = [ 1/2-1/4 0; 0 1/2-1/4];

E = kron( L, eye(n) );

A = 1/4*ones(2*n);

F = A + E;

for i = 1:n
F(i, 2*i-1) = 0;
F(2*i - 1, i) = 0;
end

Which strikes me as particularly inelegant & unnecessarily verbose (I have heard several times that if you have to use a for loop in matlab, you're probably doing it the wrong way . . .!) but I don't know how to do better.

Thanks so much for having a look.

All the best,

Rax

I think this is what you want (replace 10 with your number of rows).

for matrix:
2 -1 0 0 0 0 0 0 0
-1 2 -1 0 0 0 0 0 0
0 -1 2 0 0 0 0 0 0
...
use:

toeplitz([2,-1,zeros(1,10)])

--

Sorry didn't realize how old this post was -- must have been necro'd by someone/

Last edited:
Ah - if you were replying to my post, that's not really what I was trying to get, for N = 2, it would look like this:

0.50 0.00 0.25 0.25
0.00 0.50 0.25 0.25
0.25 0.25 0.50 0.00
0.25 0.25 0.00 0.50

Thank you for having a look!

Cheers,

Rax

Yeah, I was originally responding to the first, then saw yours, then realized you were asking a different question.

Can't think off hand of something super elegant for you.
This is best I can do off-hand, sorry :(

blocks=blkdiag(A,B);
background1=k*ones(10,10);
background2=blkdiag(-k*ones(5,5),-k*ones(5,5))
final=blocks+background1+background2;

What would be the easiest way to generalize this for arbitrary matrix size? Specifically, how can I use the blkdiag() function without having to manually insert -k*ones(n/2, n/2) over and over again?

This is what I was getting at when I was trying to define a matrix of arrays (i.e. a tensor) instead of just a matrix - but I can't figure out how to do this...

What would be the easiest way to generalize this for arbitrary matrix size? Specifically, how can I use the blkdiag() function without having to manually insert -k*ones(n/2, n/2) over and over again?
Try playing with something like...
A=-k*ones(n,n)
kron(eye(M),A)

rax,it is the same as the block diagonal. So for your last example,

B = 0.5*eye(4) + 0.25*(kron([0 1;1 0],ones(2))

will do. Basically it is the application of the same idea, you build from the smallest structure to the bigger one. Of course, there are shorter ways to do it. But it would be an overkill for these sizes.

But it would be an overkill for these sizes.
I'm just stepping out to write an exam (joy), so haven't had a chance to run this yet (hard to visualize before I do!) - but the size will grow substantially N will be on the order 100s-1000s (I hope not larger) not including the factor of 2 . . . would this change your approach at all?

Thanks very much for this, both of you - it's really helpful in terms of understanding how Matlab works. Don't know how you figured it out yourself!

Cheers,

Rax

Yes, I guess it would change the approach but I don't think that it would be something that slight modifications can not fix. When you are done, give us the actual 100-1000 matrix structure and we can figure it out together. (I am not sure what you want on the sup and super diagonals)

Thanks very much for this, both of you - it's really helpful in terms of understanding how Matlab works. Don't know how you figured it out yourself!
I've found the Octave reference manual pretty useful. It doesn't give many examples or how-to's, but you can flip to a section and get a pretty good idea what tools are available.

So the code is running alright for the time being, but apparently draws on too much memory for > 2000 x 2000 (ish).

I'm working on an in-between step now, not requires me to call on this code, run it for several values of one of the variables, generate a graph with title, axes, grid etc and save it to a specific file.

I've figured out how to do everything except the last step (viz. saving the generated figure to a *specific* file that is not simply the matlab directory). Would I simply include a change director "cd - " command *in* the .m file to do this?

Also, is there a more efficient way of creating a vector of the outputs of my function than running a for loop? i.e. is it possible to write a .m file that takes a range of values, with a specified increment and will calculate the output for each of those values *without* using a for loop?

Again (&again!) thanks for the help. I've spent so much time searching around on line to no avail.

Cheers,

Rax

Well if you are not doing complicated stuff, then, you can vectorize your variables, and send it to your function as vectors instead of computing one by one. it really depends on what you are trying to do.