This content has been marked as final. Show 21 replies
what's the link to calculating the jacobian of a real symmetric matrix?
As best I can figure out, Jacobi rotations are used to find the eigen vectors and eigen values of a symmetric matrix.
First I found this: http://en.wikipedia.org/wiki/Jacobi_eigenvalue_algorithm
But the pseudo code there kinda doesn't make sense to me and I think they have left out some steps. So I kept searching and found this:
I think that is FORTRAN and I don't really know that. Most of it makes sense and is well explained, but I don't really know what this call means:
jrotate(a(1:ip-1, ip), a(1:ip-1, iq))
I see the subroutine, further on, but I don't understand the FORTRAN notation.
a is my matrix
ip is a row in that matrix
iq is a column in that matrix.
But I don't get the part before the colon.
Or alternately if I could understand what:
wk1(:) = a1(:)
a1(:) = a1(:) - s*(a2(:) + a1(:) * tau)
a2(:) = a2(:) + s*(wk1(:) - a2(:) * tau)
means I would be set.
i used to program in fortran. but the 2nd reference i cannot understand.
your first reference i understand.
do you have questions about the first reference that might make it easier to implement?
Thanks for any help you can give me. Let's start with the "pivot" as described near the end of the Definition section. I'm taking it that this isn't the regular matrix "pivot" (diagonal). So the largest off-diagonal it is. Now in the Cost section they talk about introduction an array M (i=1…i=n-1). I think this is what is in the code as ind[k], but in their sample it goes to n, not n-1. I do get that largest means greatest magnitude, and will just keep using largest.
If M is supposed to be n-1, then I'm guessing it isn't just the largest off-diagnal, but rather just the largest in the upper triangle? But if that is so, then the code given for maxind(), doesn't really seem to make sense. The discussion page has somebody who said that code didn't work, but I'm just really stymied there. So once this is correct, I'm thinking ind would hold the potential pivot points for each row of the upper triangle.
Then there is the whole issue of this pseudo code has 1 based arrays/matrixes, whereas Flash is zero based. So I keep getting confused on which loops need to be adjusted. So this is also makes it hard for me to understand the first part inside the while statement. It looks like it is going through all the potential pivot points to find the largest one. But the way the pseudo code is written it seems to suggest that you can't pick one is row 0 or row n-2 or row n-1?
Anyways, I've attached the code I've come up with so far. I renamed a few of the variables just to help me follow it. Also instead of a while I replaced it with a for loop. (Since it wasn't converging I got tired of waiting for the slow script dialog!) If there is anything you can suggest that would be great.
Also the text says that the upper triangle will be destroyed, but that the original Matrix can be restored by copying the lower triangle. But isn't the diagonal also going to be changed?
i'll take a look at your code, but that's going to take a while to go through.
just a quick note: he's trying to perform various rotations to get the part of the matrix above the diagonal to all zeros. that would convert it to a triangular matrix and because the eigenvalues of a triangular matrix are its diagonal elements, the problem (finding the eigenvalues) is solved.
that is guaranteed to be doable because the matrix is symmetric. and yes, the diagnoal elements will be changed with his rotations. if they weren't changed it would be trivial to list the eigenvalues of a symmetric matrix.
That is more or less what I thought, so that is why I'm confused by some of the statements made in the narrative and then the corresponding (or actually not corresponding) instructions in the code.
I spent several hours last night and several more today making my head hurt a bunch with this. Any help or insights you can provide would be very helpful.
i checked your code and i think there are some errors. however, correcting them doesn't yield the correct answer for the matrix listed on wikipedia nor for a 3x3 matrix that i tried.
so, i also conclude the algorigthm is incorrect. i'm checking now to program an accurate representation of jacobi's algorithm.
i think i just wasted 2 hours of my life. and i don't even want to check this any further in case i find a problem. i checked it against 3 symmetric matrices using mathematica and it was accurate.
Wow. Thank you so much. You are awesome.
I hate to look a gift horse in the mouth, but…this returns the eigenvalues, but I'm also trying to get the eigenvectors. I'm certainly willing to give it a go on my own, but do you have any hints as to where I should look?
All the examples I've found are good for if I'm writing it out, but they all involve elementary row operations which I don't understand how to program.
it never hurts to ask. and i had a good night's sleep anyway.
plus, there's no extra work to do. v is the eigenvector matrix and i'm pretty sure the eigenvectors are arranged as columns in v.
p.s. i didn't use tol which should be the tolerence or accuracy of calculation. like your attemps, i didn't want to deal convergence problems because of more basic formula issues.
You are so amazing on so many levels. Now that I've had a good nights sleep too I see that is in there. Dealing with the tolerance is something I can do.
Wow. Wow! and wow.
Do you have any advice for me on the other part of the problem? Just the general idea of making a Matrix class? That I can do, just if there are anythings you would suggest I look at or any ways that I should think about structuring it?
Thank you very much!
i have to run to my office. will be back in an hour or two. i'll post my ideas (if i have any worthwhile) when i return.
here's the code for using the tolerance and returning the eigenvectors:
Thanks again, now the questions begin! :) These shouldn't be too hard.
If the answers on the wikipedia page are to be assumed correct, it seems that some of the eigenvectors come out with the wrong sign. I'm thinking that may have something to do with Flash's upside down coordinates? Or perhaps how the Math.atan2() functions. I can keep poking around, but if you have any suggestions.
The next thing is that yes the eigen vectors are in V, but they go across the rows. I would like them to go down the columns. I can easily do a transpose before returning them, but if there is a quick "swap the i,j,k array indexes and your good as gold" that would be great.
And wow this is so fast. Including the tolerance cuts it down to just 4 to 12 iterations for 5 to 150 or so milliseconds for the 4x4 and 9x9 matrices I've tried. And that even includes where I've added copying the original matrix so that it won't be destroyed.
i don't see any problem with the eigenvectors. the eigenvectors found by the script agree with those found by mathematica.
mathematica does appear to standardize (but not normalize) them so i added a matching standardization routine to the code to make it easy to compare the a.s. results with mathematicas.
the speed of this code is because jacobi's method converges quickly not because of any programming skills on my part.
below is the standardization code and transpose code (which you should have in your matrix class):
actually, you should adjust the transpose code so it doesn't replace the original matrix.
Thanks again. If it agrees with Mathematica that is good enough for me!
I'll incorporate the other changes too.
do you have any application in mind for your matrix class?
I'm trying to see if Flash can handle facial recognition and identification. Mostly just as a lark.
I'll let you know if I'm able to get it to work.
whoa, i wouldn't have suspected that. are you using a webcam and the bitmapdata class to identify features?
Yup, that is more or less it. Once I've got the basics down I have a few things I want to try with it. I'll keep you posted.