Tuesday, January 24, 2012

Symmetric Matrix

I had a 4x4 symmetric matrix to be entered in Octave and I was lazy to feed redundant elements :), convincing myself that it is waste of time. Well, I spent more than that 'to be wasted' time browsing for how to enter symmetric matrices economically in Octave. Believe me, it was worth it. You never get bored fiddling with matrices and of course, Octave/Matlab.

What I found was not a single step solution, for that matter I could not get a single step way to do it at all. I was impressed by this solution for its elegance. Not just that, this solution caters to my requirement of feeding unique elements (either upper triangular or lower) exactly once. Okay, now on to the steps:

1) Say your matrix is as below:
    16    4    8    4
    4   10    8    4
    8    8   12   10
    4    4   10   12

Obviously symmetric - let's call this A


2) V=[16 4 10 8 8 12 4 4 10 12];

V holds all the upper (lower) triangular elements one column (row) at a time.

3) Create a triangular matrix of 1 of the same size as your A matrix as follows:

M = triu(ones(4))
M =

   1   1   1   1
   0   1   1   1
   0   0   1   1
   0   0   0   1

4) Replace the ones with elements from V as follows:

M(M==1)=V
M =

   16    4    8    4
    0   10    8    4
    0    0   12   10
    0    0    0   12


The M==1 test returns the indices of all places that has 1 in it, and the assignment takes care of replacing them with elements from V.

5) Add M with itself but after transposing and taking only the lower triangular portion (note M is an upper triangular one):

M = M + tril(M',-1)
M =

   16    4    8    4
    4   10    8    4
    8    8   12   10
    4    4   10   12

As you can see, M is A.

-1 in the tril command makes sure that you are leaving the main diagonal when retrieving the lower triangular portion of the transposed M.


So, as it looks it is a 4 step algorithm, but if you can create a script with these steps in it, it is only one step.

Now to credits: I got this tip from one of the threads of mathworks.com :), where else?!

No comments:

Post a Comment