Getting the corresponding eigenvectors in arpack

General Usage
1

I’m using the arpack eigs function to retrieve the 3 largest eigenvalues in an 8x8 matrix ,using the following call

eigs(a, nev=3)

works great and i get a vector of 3 eigenvalues.

However i get an 8x3 matrix of eigenvectors.

This question probably does not make sense (go easy on me it’s been a while since i’ve had to work eigenvalue/vector problems and never a “reduced” problem), how do i select the 3 eigenvectors that correspond to the eigenvalues i’ve chosen ?

edit: ouch. I need to rethink this. The paper i’m working from is using a Lanczos iteration to generate the reduced matrix “T”. it then applies an eigendecomposition to the matrix T which provides the eigenvalues and eigenvectors of the reduced system.

this seems odd since i thought that the whole point of Lanczos was to provide the “dominant” eigenvalues of a system directly.

anyway, if you have any experience in this sort of problem feel free to point me in the right direction !

thank you

2

For these very small matrices you can simply solve the dense eigenproblem, i.e. you get the full decomposition with all (here 8) eigenvalues and vectors. You can then sort your eigenvalues using sortperm and access the largest eigenvalues and eigenvectors with the corresponding indices. The eigenvectors, are stored as matrices with the second dimension being the index of the eigensolution/ corresponding eigenvalue. This is identical to the solutions of eigs in Arpack/ArnoldiMethod.jl/etc…

julia> using LinearAlgebra

julia> a = rand(8,8);

julia> evals, evecs = eigen(a);

julia> evals
8-element Vector{ComplexF64}:
   -0.804790799527211 + 0.0im
  -0.6140770220738357 + 0.0im
 -0.20256099079070283 - 0.7372771202355107im
 -0.20256099079070283 + 0.7372771202355107im
 -0.07136837472348928 + 0.0im
  0.12426312411688796 + 0.0im
   0.6786223201706776 + 0.0im
   4.2558694471181875 + 0.0im

julia> per = sortperm(abs.(evals),rev=true);

julia> evals[per[1:3]]
3-element Vector{ComplexF64}:
   4.2558694471181875 + 0.0im
   -0.804790799527211 + 0.0im
 -0.20256099079070283 - 0.7372771202355107im

julia> evecs[:,per[1:3]]
8×3 Matrix{ComplexF64}:
  -0.32752+0.0im    0.524908+0.0im  -0.0116975-0.465659im
 -0.404005+0.0im  -0.0655954+0.0im  -0.0489684+0.442762im
 -0.287834+0.0im   -0.352892+0.0im    0.454449+0.141423im
 -0.400029+0.0im   0.0318565+0.0im   -0.142492+0.013602im
 -0.236407+0.0im    0.547071+0.0im   -0.154853-0.0438481im
 -0.306862+0.0im   -0.201882+0.0im   -0.478091-0.0im
 -0.380962+0.0im   -0.479747+0.0im  -0.0821356+0.00777412im
 -0.437556+0.0im   -0.156271+0.0im    0.276123+0.00721887im

For your question on the Lanczos iteration, I believe that this is just the underlying algorithm of the eigensolver the paper you refer to has used. When the matrix a is symmetric, the Lanczos iteration is the preferred algorithm to solve for eigenvalues/vectors. You can let julia know that your matrix is symmetric by calling a=Symmetric(a) before you call eigen or eigs.

Does this help?

3

?? The eigenvectors are the 3 columns of your 8x3 eigenvector matrix, returned by eigs.

Oh, I think I see. You started with a big matrix A, then reduced it to an 8x8 matrix T using Lanczos (implemented by you, not by Arpack), and found some eigenvectors y_k and eigenvalues \mu_k of T. Now you want to find the corresponding (approximate) eigenvectors of A?

The corresponding (approximate) eigenvectors of A are x_k = Q y_k, where Q is the orthonormal basis constructed during the Lanczos procedure. (Read up on the Rayleigh–Ritz method.) However, the simplest implementations of Lanczos don’t store Q — they only use a 3-term recurrence to construct T (storing only the most recent 2 columns of Q), which means you can’t (easily) get the eigenvectors anymore. You have to modify your algorithm to save all of the previous q vectors to keep track of the orthonormal basis explicitly, more like Arnoldi.

However, there’s usually absolutely no reason to implement Lanczos yourself except as a learning exercise. (Especially since it’s easy to get “ghost eigenvalues” if you implement it naively.) Just use the one in Arpack.jl (passing Hermitian(A)), or in KrylovKit.jl (pass ishermitian=true), or use LOBPCG in IterativeSolvers.jl (a different iterative algorithm for Hermitian problems).

4

A is definitely not symmetric.

I took the idea that you proposed, @fgerick, and solved the dense problem taking the 3 dominant eigenvalues and vectors. It gives a reasonable result but i don’t get a good approximation. After some thought i realized that’s probably the correct behavior and this approach won’t work in this particular case, i.e. the end goal is to construct a Pade polynomial which represents the transfer function of the system.

I think it’s equivalent to taking an 8th order polynomial , chopping down to a 3rd order and expecting it to be an accurate approximation - it won’t be.

The big matrix A is the 8x8 matrix. Instead of solving for all the eigenvalues/vectors to construct the frequency response of the network, it is reduced to a 3x3 matrix using Lanczos, i.e. the resulting T matrix is a 3x3.

I have implemented the (naive) Lanczos algorithm AND have used Arpack. I was not clear about that.

The toy lanczos process works, i.e. it gives me a reasonable T matrix. I then run eigen on the T matrix and i get eigenvalues/vectors that work, but not well. Which then made me think that my toy lanczos does NOT work, so i went to Arpack/KrylovKit.

When i use Arpack i get the 3 eigenvalues, which correspond very closely to the values i get from the toy Lanczos, but i’m unclear on which eigenvectors are associated with the eigenvalues, because I have an 8x3 matrix but only need a 3x3 result (the eigenvectors are needed to calculate the residues of the polynomial).

I should point out that I’ve used a different method (AWE) to calculate the polynomial equivalent and a 3rd order (Pade) polynomial can be made which is a very accurate reproduction of the original system. So i’m not chasing a red herring.

I’m most definitely in over my head so i’m camped out in chapter 9 of Golub and van Sloan trying to figure it out :slight_smile:

I just looked at the Rayleight Ritz method and that looks more directly applicable to what i’m trying to do rather than using Lanczos. The papers i’m referring to are quite old (90s) and maybe using Lanczos to do this is not the best approach anymore.

I’m not sure why you suggest using Arpack and passing in Hermitian(A) instead of just A.

Thanks very much for your help !

5

This makes no sense — if it’s not Hermitian, you can’t use the Lanczos algorithm [*], and you need to use a non-symmetric algorithm like Arnoldi. (Lanczos is the special case of Arnoldi for Hermitian operators.)

[*] Unless A is Hermitian under some modified inner product \langle x, y \rangle = x^* B y, i.e. if B A is Hermitian for some Hermitian positive-definite B.

Lanczos (or Arnoldi) is how you get the subspace basis that you plug into Rayleigh–Ritz.

Because you said you were using the Lanczos algorithm, which implies that your A is Hermitian, in which case it is beneficial to tell Arpack that your matrix is Hermitian by wrapping it in the Hermitian type.

6

The matrix is not symmetric.

aha. There is a paper referenced which contains the method the authors say they actually used - “An implementation of the look-ahead Lanczos algorithm for non-Hermitian matrices.”

now i understand my confusion. They sketched out the naive Lanczos[*] algorithm in the paper which it turns out is a really bad idea if someone like me is reading it because the standard naive method doesn’t even apply to a non-symmetric matrix (I just got to section 9.4). I implemented it and it gave me reasonable but not very good results (the fact that it didn’t just fall flat on its face seems surprising).

I need to dig into the paper trail some more and see what i can figure out. I’m quite sure i’ve seen some other papers on this subject which explicitly talk about using the Arnoldi method.

My goal is definitely to use Arpack or KrylovKit, but in the current situation I don’t see how I can because of my confusion around the eigenvectors. perhaps i can figure it out from some other papers.

Thanks very much again.

[*] what i thought was the plain Lanczos algorithm appears to be from the paper I referenced, i.e. it is the unsymmetric implementation. They referred to the original Lanczos paper and the original paper also appears to contain an unsymmetric algorithm. This is going to take a bit to unwind…

7

Oh right, there is indeed a non-Hermitian variant of Lanczos. It’s not very commonly used AFAIK, in contrast to Arnoldi. (Maybe because it sounds a little less reliable than Arnoldi, because it can break down, added to the usual Lanczos problem with ghost eigenvalues.) Most references to Lanczos that I’ve seen are referring to the Hermitian version.

If you have an 8x8 matrix, you should not be using an iterative solver like Arpack or KrylovKit. You should just call eigen to get all 8 eigenvalues and eigenvectors. It will be way faster.

Once you have the eigenvectors you can project the matrix on to whatever subspace of them you want. It’s not clear to me what you actually want mathematically at the end, though.

The 3x3 “T matrix” that you get from a Lanczos procedure is similar (via a change of basis) to just the corresponding portion of the diagonalization of A in the eigenvector basis, so anything you got from the Lanczos T matrix you should also be able to get directly from the eigenvalues and eigenvectors.

1 Like
8

What i want at the end is the poles and residues of my system.

The goal of this methodology is precisely to avoid calculating all of the eigenvalues and eigenvectors. The idea is that i have a circuit matrix that is quite large but the frequency response of the circuit is such that only the first 8, 10, 12 etc… poles really matter, therefore you don’t need to call eigen on the large matrix.

Obviously my 8x8 matrix is a toy problem used by me to make sure it all works, a realistic problem would very easily be 100x100. And yes practically speaking, calling eigen even on 100x100 matrix still probably isn’t that much of a problem. However, I have a particular application in mind where i do need to calculate the minimum polynomial that will give me a good approximation to the full solution.

Your comment about projecting the matrix onto my desired subspace in really interesting to me. The problems that i’m currently trying to solve, even the large problems, are still sized such that running eigen directly would not be an issue. What is an issue is having a resulting polynomial of degree 50 or 100.

So if i had a method to take the full set of eigenvalues/eigenvectorrs and create my reduced result i would do that.

9

For what it’s worth, I have posted an implementation of the Look-ahead Lanczos procedure here: https://github.com/JuliaLinearAlgebra/IterativeSolvers.jl/pull/321

Also note that you don’t necessarily need look-ahead to make it work for an unsymmetric matrix: https://epubs.siam.org/doi/abs/10.1137/0719031?journalCode=sjnaam . This is the variant that you can find in the guts of QMR iterative solvers in KrylovKit, etc.

1 Like
10

I still don’t fully follow the problem. You could just calculate the dense eigen solution for the 100x100 matrix (that takes a few milliseconds on my laptop), then you select your largest eigenvalues/vectors. And if you then want to take only a few coefficients from the eigenvector (I assume they’re coefficients of monomials/polynomials!?), you just take the x number of largest components of the eigenvector? That’s not really a mathematically clean method, but it seems simple to me.

11

The end result of the process is calculating the poles and residues of the polynomial, H(s).

if you use the full solution, which i have done , you get a “n” degree polynomial. n is not necessarily the size of A, which has to do with how the circuit equations are used to generate the matrix. In my example n =6, not 8.

However the goal of the problem is to create a lower degree, “q”, polynomial which is the equivalent of the full solution, in this case i am using 3.

Simply keeping the largest 3 eigenvalues (the poles) and manipulating the eigenvectors to provide the values i need for the residues gives me a reasonable response but one that is not very accurate.

Admittedly i cannot explain this mathematically because i’m still trying to catch up on exactly how this works, but intuitively i think it’s because you are simply truncating the polynomial and expecting that to provide an accurate response. Clearly you would need to recalculate the coefficients of the polynomial of lower degree to have any hope of it matching the higher degrees polynomial response.

12

Is this similar what you are trying to do: Model reduction methods based on Krylov subspaces? Your original post emphasized getting eigenvalues and eigenvectors. Computing eigenvalues is by far the most common use of Lanczos, and “Lanczos” unqualified is usually going to be interpreted as symmetric Lanczos. For model reduction, eigenvalues aren’t the central point, but instead you use non-symmetric Lanczos as a tool to compute Padé approximants of a transfer function.

The above survey article is quite general and presents the block form of the algorithm, which is appropriate if you have a MIMO system. If it is SISO, it’s a bit simpler and a standard nonsymmetric Lanczos is what you want, possibly with look-ahead for stability, although you might get lucky with your specific problem and not need it. I know about these algorithms in a shallow way, but I’ve never really used them and don’t have any practical experience. And I’m not sure what you want to do with this after nonsymmetric Lanczos. You do have an expression for the transfer function at that point and the poles are indeed the eigenvalues of T. You shouldn’t need the eigenvalues/eigenvectors of the original matrix.

13

Yes, that’s what I’m trying to do. Thanks for the link I hadn’t seen that paper yet.

My system is SISO.

if you do a QR on the T matrix you can simply use the resulting matrices to calculate the transfer function without calculating the poles/residues.

I haven’t tried that yet. I should probably do it to see if it provides any insight as to why it’s not working.

I forgot to mention. AWE works really well as long as long as the order of the polynomial is not too large. “too large” is problem dependent, but is generally 10-15.

The idea behind the Krylov subspace methods is that it allows generating much higher degree polynomials.

14

I’m not sure I understand what’s not working for you. ARPACK is going to attempt to give you eigenvalues and eigenvectors of the large matrix. That was your problem in your original post with the 8\times 3 matrix of eigenvectors. You’re perfectly right that you can use QR factorization to evaluate the transfer function quickly once you have T. So if you have done the nonsymmetric Lanczos and obtained the correct T for whatever order you determine gives an acceptable reduction, you should be set. If there is some more subtle problem with the T you computed, it’s not clear to me how you are deciding that it worked or didn’t work.

15

I’m currently deciding that it did not work by evaluating the resulting polynomial, which does not give an accurate representation of the system.

One thing that i can check is, does my Lanczos give me the same eigenvalues as Arpac ? The answer is yes for the first two and no for the 3rd.

Therefore it appears I’m NOT constructing T correctly, which is hardly surprising since it’s non trivial to implement and i’m using the “simplified” version found in the original paper (although it does appear to be a simplified implementation of the unsymmetric version). It’s essentially simplified because it does not implement lookahead, which should be fine because i can clearly see the “delta” term being calculated is not anywhere near zero.

Humorously this entire discussion started precisely because i was using Arpack to try and check my Lanczos results and got stuck because i didn’t know how to use the 8x3 matrix of eigenvectors. I have 3 vectors of length 8. You would think that the first 3 elements of each column would be what I want and I would be done. However if I do that I don’t get a correct polynomial either.

I’m now beginning to think that the way i’m using the results of Arpack are correct and the real problem is that i’m missing some subtlety in the paper as to how i’m calculating my polynomial values. Hence the reason I was going to pursue the QR method, as an additional check.

16

Not even close.

17

ok. Reading back through your comments i see that your plan would be to simply take the full eigenvector/eigenvalue results and reduce them to my q=3 case.

Any hints as to where i might read up to understand how best to do that ?

18

This may be of interest to future travelers of this thread…

https://www.math.ucdavis.edu/~freund/BANDITS/

19

I’m wondering about “evaluating the resulting polynomial.” The algorithm is for approximating a rational transfer function, so I’m not sure about that polynomial.

ARPACK is not going to be a useful check. The model reduction algorithm uses Lanczos to compute a reduced model with a transfer function that is a Padé approximant to the original transfer function. Lanczos probably will generate some good eigenvalue approximations in T along the way to computing that approximant, because that’s something it does, but that’s not a necessary part of computing the Padé approximant and there is no particular reason why you need all the eigenvalues of a smaller T to match eigenvalues of A or eigenvalues computed using ARPACK. You are using Lanczos for reasons that are not so directly tied to computing or approximating eigenvalues.

If you do want to check that the Lanczos procedure is working, you can use it to compute the full 8\times 8 tridiagonal and then compare those eigenvalues to those of A. Those have a decent chance of matching reasonably well on a small problem, with any differences being due to numerical error that with luck won’t be too large on your 8\times 8 problem. Or you could check the underlying relations that are supposed to hold in the nonsymmetric Lanczos algorithm (e.g. biorthogonality). You could check that the reduced transfer function is a good approximation to the original around the point you are expanding at. Or you could check that you in fact go the Padé approximant at the end, which is what is supposed to be achieved. There are a lot of things you can look at. But comparing eigenvalues against ARPACK won’t work.

20

i was being sloppy, it is a rational polynomial. To be more specific

\Sigma_{j=1}^{q}\frac{k_j}{\sigma - p_j}

The p_j are exactly the eigenvalues. The k_j are calculated using the eigenvectors, hence my interest in the eigenvectors.