3

I'm writing an algorithm to find the eigenvalues and eigenvectors of a positive definite matrix with the power iteration method. I know it's very crude, of course, and there are better methods, but this is just a trivial application and I don't want to go much beyond (nor do I have access to any libraries, the language doesn't have any). I'm finding the eigenvalues with a power iteration, then at the end of each I redefine the matrix as:

$ \mathbf{M}^{(n)} = \mathbf{M}^{(n-1)}-\lambda_{n-1}\mathbf{e}_{n-1}\mathbf{e}_{n-1}^T $

to remove the eigenvalue, and repeat the process. I get pretty excellent values of the eigenvalues - they match with the solution I can get with NumPy up to 1E-6 precision, easily. However the eigenvectors for all but the first one or two eigenvalues are a complete mess. I perform a Gram-Schmidt orthogonalisation on them after the power iteration finishes, and I even check that they return the correct eigenvalues with the original matrix as their Rayleigh quotient - they do - but still, they're very different from the ones I get from NumPy. What could I look into? Is it just a matter of numerical noise, and there is no chance to improve unless I move to better algorithms, or am I missing something obvious?

Rodrigo de Azevedo
  • 23,223
Okarin
  • 153
  • 7
  • one thing needs confirmation. You say positive definite, do you mean real symmetric (with positive eigenvalues) – Will Jagy Nov 27 '19 at 23:20
  • No, it's a matrix of a system of coupled oscillators. I know it's positive definite because the eigenvalues are all $\omega^2 > 0$, and it's all real, but it isn't symmetric. – Okarin Nov 27 '19 at 23:28
  • Alright, maybe someone will be able to help. Meanwhile, I found out who first wrote down that switching Hamlet and Othello into each other's story would result in both winning easily; about 1905, Andrew Cecil Bradley. There is a slightly longer mention in lectures by Harold Bloom about 1988. https://www.britannica.com/biography/A-C-Bradley – Will Jagy Nov 27 '19 at 23:40
  • That's... huh... interesting.

    Anyway I suspect it's really just numerical noise. I tried putting the eigenvector for the smallest eigenvalue from my program into NumPy, the resulting Rayleigh coefficient is almost identical, but just a bit bigger. However my program misses that and finds that it's perfectly equal to the eigenvalue. It might be because it doesn't have double precision (my program is written in Godot Engine, a game engine with its own scripting language, not exactly designed for scientific computation).

     – Okarin Nov 27 '19 at 23:40
  • I use C++ for most things, it is one of just two computer languages in which I was able to take a course. The person I knew who was doing scientific computation full time really loved Java. – Will Jagy Nov 28 '19 at 01:24
  • By the way, usually if you say a matrix is positive it is automatically assumed it is also hermitian. – lcv Nov 28 '19 at 07:46
  • @WillJagy yes, obviously if this was a serious application I'd be using C++ or even better, Python... but my goal here was exactly to make an algorithm that could diagonalise matrices in this specific language, so I've got to work from scratch and with some really annoying constraints. – Okarin Nov 28 '19 at 10:13
  • @Icv however I'd say that's a stretch... positive definite means all the eigenvalues are real and positive. It doesn't need to be Hermitian or symmetric for that to be true. – Okarin Nov 28 '19 at 10:15

2 Answers2

2

The missing point is that the matrix is not symmetric (usually implicitly assumed with positive definiteness). If it is not symmetric, then the eigenvectors are very likely not orthogonal so there is no reason to assume that $M^{(n)}$ and $M^{(n-1)}$ share the same eigenvectors.

What you need to do is to use the update $$ M^{(n)}=M^{(n-1)}-\lambda_{n-1}e_{n-1}f_{n-1}^T, $$ where $f_{n-1}$ is the left eigenvector corresponding to the eigenvalue $\lambda_{n-1}$ normalized such that $f_{n-1}^Te_{n-1}=1$. It can be computed by the same power method applied on the transpose of $M^{(n-1)}$.

Algebraic Pavel
  • 23,818
  • Thanks! Yes, I just realised that was my mistake too. The algorithm works for symmetric matrices, so that's definitely it. – Okarin Nov 28 '19 at 10:10
  • 1
    You're welcome. I'm still a bit puzzled why the eigenvalues are correct with your original process but maybe it just works for some reason :-) – Algebraic Pavel Nov 28 '19 at 10:16
  • I found out. What I've been using is called the Wielandt deflation, and it works for a number of different choices for the vector on the right; it doesn't even need to be an eigenvector at all! In fact my original choice is considered one of the optimal ones for the eigenvalues; it also preserves the left eigenvectors, but not the right ones, which was my problem. – Okarin Nov 28 '19 at 11:23
  • @Okarin Ah right, thanks for the info. – Algebraic Pavel Nov 28 '19 at 12:31
1

There are two problems at play here.

First: In a power iteration method you usually don't redefine your matrix by getting rid of the dyadic problem of the eigenvectors. It's more common to simply subtract the projection to already found EVs from your current iteration. You do $\tilde{v}_i = v_i - \lambda_1 e_1 e_1^Tv_i$ and continue with $\tilde{v}_i$. This is mathematically equivalent to your approach but does not rely on the calculation of dyadic products and does not destroy any sparsity you have in $M$.

Second: The power method works by projecting your current iterative into the EV-decomposition $$v_i = \sum \limits c_k e_k$$ where $e_k$ are the EVs. One iteration now gives you $$ v_{i+1} = \sum \lambda_k c_k e_k $$ i.e. every coefficient is multiplied with the eigenvalue. Over time, the largest eigenvalue will dominate and be the only remaining one, numerically speaking. Since you have a matrix with complex oscillators, at one point you don't have a single largest Eigenvalue anymore. At this point, the algorithm can do a lot of things, including convergence to any of them or periodically switching from one to the other. Since you are going to find complex EVs, using a real $v_i$ will never give you a good result for an EV.

Edit: OP clarified, that there are no complex eigenvalues. Theory about values with the same magnitude still holds.

Laray
  • 2,398
  • The sparsity isn't a major worry because the matrix is small - as I said, this is basically a toy application, nothing serious.

    About the rest, as I said, the eigenvalues are real, positive, and I find them all to excellent precision. It's the eigenvectors that are wonky, and I think I realised why: I'm too used to working with symmetric matrices and mixed up left and right eigenvectors. When I try the program with a symmetric matrix, it works well enough.

     – Okarin Nov 28 '19 at 10:03
  • If you don't use a symmetric matrix, you can not guarantee real eigenvectors. Secondly, you talked about coupled oscilators, that I always read as "has complex Eigenvalues" and corresponding eigenvectors. – Laray Nov 28 '19 at 10:42
  • These are coupled, non-damped, non-forced oscillators. I don't see why they should have complex eigenvalues? The eigenvalues are literally just the $\omega^2$. In practice I know for sure the matrices I tested it on were positive definite and had real eigenvectors. – Okarin Nov 28 '19 at 11:21
  • Oscillations occur in control-technology, when the system matrix has complex eigenvalues. That's where my association comes from. If that's not your case, I gladly admit, that I was wrong... – Laray Nov 28 '19 at 11:32
  • I think for the kind of system I'm describing that would be unphysical. Control technology usually involves PID, which means you got effectively both first and second time derivatives... that's equivalent to a damped oscillator, and yes, those ones have complex eigenvalues (the imaginary part corresponds to the decay). – Okarin Nov 28 '19 at 11:42