[Trilinos-Users] PyTrilinos + Anasazi - Generalized eigenproblem

[Marko Budisic]

Dear Chris,

thank you for your prompt response. Let me state at the beginning that I 
will usually be working with N on the order of 1e3 up to possibly 1e5. 
So I thought Anasazi would be a good pick. I was using a 5-dim problem 
because I could easily (and quickly) generate something that resembled 
my intended application and check whether I was using the code properly. 
Now I understand that some of the problems were rooted in the small size 
of the operator. I'm including the rest of the responses where 
appropriate below.

Chris Baker wrote:
> Dear Marko,
> 
>     1) It seems that the setM() operation is ignored by the code, i.e.
>     setting setM(D) still computes just the regular eigenproblem, instead of
>      the generalized one.
> 
> 
> This should not be the case. I'll look into it some more. To be clear, 
> you are saying that you get the exact same eigenvectors and eigenvalues 
> whether or not you include the call to setM() ? My recollection of the 
> code is that it should have the intended effect, of scaling the 
> eigenvalues (the resulting eigenvectors will probably not be unit 
> length; I don't recall).

Correct. Even after passing a non-identity matrix to setM(), I get an 
eigenvector and eigenvalue of the non-generalized problem (due to 
problems below I couldn't compute any other eigenvectors, I will try 
generating something larger and let you know). Since I was using the LM 
and LR settings to check this, I get the evector corresponding to evalue 
4.3... The eigenvector is of unit length.

> 
>  
> 
>     2) Whenever I try to set Nev > 1, e.g. Nev = 2, the solver complains
>     about orthogonality constraints producing the following error:
> 
>     "Anasazi::BlockKrylovSchurSolMgr: Potentially impossible orthogonality
>     requests. Reduce basis size."
> 
>     I am guessing there are some issues with my setting the blocksize. I am
>     passing a randomized storage multivector to setInitVec(). I have tried
>     setting the size of multivector to both 1 and to Nev (the number of
>     eigenvalues requested) and in both cases I get the above complaint.
> 
>     Could anyone shed some light on how one would go about figuring the
>     right blocksize for this problem and in general?
> 
> 
> The operation of this solver is that blocks of vectors (of size "Block 
> Size") are added to an orthonormal basis up to a certain basis size 
> ("Block Size" times "Number of Blocks"). At that point, the basis is 
> restarted. However, the restart must be large enough to capture at least 
> the number of desired eigenvalues ("NEV"). Lastly, fundamental linear 
> algebra dictates that the size of the basis cannot exceed the length of 
> the vector space (for this problem, 5).
> 
> The default number of blocks for this solver manager is 3. This is the 
> minimum that you can use for this solver. Typically, one uses a 
> significantly larger basis. In your case, the block size and basis size 
> depend on how many eigenvectors you want to compute, but my advice is 
> the following: if you will only be solving 5x5 eigenvalue problems, you 
> don't want to use Anasazi. In particular, a current limitation of the 
> BlockKrylovSchur eigensolver is that the basis can only be run up to 
> size n-1, where n is the size of the linear operator (in your case, 
> n=5). This makes it impossible (even using auxilliary vectors) to 
> compute all eigenvalues of a sparse operator.

So if I had a larger operator size (say 1e3 dimensions), to compute the 
first say 10 evalues/evectors using the default options, it should be 
sufficient to pass 10 to setNev() and pass a multivector of size 1 to 
setInitVec()?


>  
> 
>     3) I am setting only "Block Size", "Which" and "Verbosity" options for
>     the parameter list passed to BlockKrylovSchurSolMgr. Should other
>     options relating to blocksize be set too? Should I set even the "Block
>     Size" manually?
> 
> 
> You are welcome to manually set the block size, but for this problem, 
> you can't. Block size larger than one will build a basis larger than 6, 
> which is not possible for a 5-dimensional vector space.
>  
> 
>     4) Should AuxVec also be set prior to solving the problem/is there a
>     benefit for doing so? If yes, how would one go about choosing the size
>     of this multivector and should it be initialized in a certain way?
> 
> 
> Auxilliary vectors define a space which the solver should avoid. This is 
> useful for restarting the eigensolve if some number of eigenvector are 
> already known. I apologize; this is not documented very well.

I have just remembered that I tried passing a randomized MultiVector of 
size 1 to setAuxVec() just to see what happens for the above problem. 
The code seemed to freeze and after 10 seconds or so I killed the 
process. I understand that this might again be due to the small size of 
the problem I used but perhaps the code could made to fail with some 
error in future, to facilitate debugging.

>  
> 
>     In the end, I am including my python code for this problem. At the very
>     end, numpy.linalg is used to verify the results obtained by the Anasazi.
> 
> 
> If you are going to continue solving eigenvalue problems on the order of 
> tens or hundreds of variables, or if you are going to compute a large 
> proportion of the eigenvalues, I suggest you use dense eigensolvers like 
> those provided by numpy.
> 
> Please respond with more questions, or if we can provide any other advice.
> 
> Best regards,
> Chris Baker
> 

Thank you for the detailed explanation of how the code works. Perhaps 
some of it could be added to the documentation in upcoming versions of 
Anasazi, I believe I understand how the code works much better now.

Best,
Marko Budisic