[Trilinos-Users] Anasazi Solver Managers Unconverge where MATLAB converges.

[Mehmet Salih YILDIRIM]

Hello everybody!

As some of you already know, we'd been trying to implement the distributed
version of normalized
cuts<http://www.cs.berkeley.edu/%7Emalik/papers/SM-ncut.pdf>image
segmentation algorithm. We decided to use anasazi to solve the
eigenproblem since our implementation is in c++ and anasazi provides
eigensolvers which may process distributed sparse matrices.

However we encountered a few problems, as I alread stated before. Even
though we did what the people who replied to our email told us, the problems
insist. This time, however, I have some extra information, and results of
some experiments:

1st: I transformed the generalized problem

(D - W) * X = \lambda * D * x
>

where D is diagonal and positive (a diagonal entry in ith row is the sum of
all entries in ith row of matrix W), into a standard problem as follows:

invSqrtD = inv(sqrt(D))
> sqrtD = sqrt(D)
>
> (invSqrtD * (D - W) * invSqrtD) * (sqrtD * X) = \lambda * (sqrtD * X)
>

so if we say

A = invSqrtD * (D - W) * invSqrtD
>
Y = (sqrtD * x)
>

it looks like

A * Y = \lambda * Y


which is a standard problem.

However, BlockDavidson, BlockKrylovSchur, LOBPCG and RTR(both IRTR and SRTR
modes) solver managers returned unconverged.

After this, I thought that i fill the matrix W with wrong values. So I
exported the matrix A (that is, invSqrtD * (D - W) * invSqrtD) using
EpetraExt::RowMatrixToMatrixMarketFile to a file, then I loaded the matrix
into MATLAB using the mmread function, and called eigs function to solve the
standard eigenproblem. MATLAB solved the eigenproblem.

So I think that I am initializing the solvers with wrong parameters. Here I
give some details:
the matrix W (as well as D, A and all others) are 3072 to 3072 (or larger)
matrices. Parameters to solvers are:

    Teuchos::ParameterList parameterList;
>     parameterList.set("Maximum Restarts", 200); //Maximum Iterations is
> also 200 in some other solvers
>     parameterList.set("Block Size", blocksize); //block size is 100
>     parameterList.set("Num Blocks", 100);
>     parameterList.set("Which", string("SR"));
>     parameterList.set("Convergence Norm", string("2"));
>     parameterList.set("Convergence Tolerance", .05); //accuracy is not
> important.
>     parameterList.set("Verbosity", Anasazi::Warnings + Anasazi::Errors +
> Anasazi::FinalSummary);
>

The final summary, for example from BlockDavidson is:

                          BlockDavidson Solver Status
>
> The solver is initialized.
> The number of iterations performed is 19899
> The block size is         2
> The number of blocks is   100
> The current basis size is 200
> The number of auxiliary vectors is 0
> The number of operations Op*x   is 79600
> The number of operations M*x    is 0
> The number of operations Prec*x is 0
>
> CURRENT EIGENVALUE ESTIMATES
>           Eigenvalue         Residual(M)         Residual(2)
>
> --------------------------------------------------------------------------------
>         5.499616e-16         not current        1.415645e-15
>         1.908911e-03         not current        1.559112e-15
>

Another interesting point is, if Convergence Tolerance is set to 0.1,
solvers perform 0 iterations and they return Anasazi::Converged.

If one will be interested, the Matrix Market file of the matrix A (for which
standard eigenproblem should be solved) can be downloaded
here<http://science.isadamlari.org/files/Amatrix.dat>.
And anyone interested, may download the full source code
here<http://science.isadamlari.org/files/distncut.tar.gz>.
Please note that the source code is not complete yet, we are still writing.

Any help would be appreciated since this is our senior project and we have
to finish it in a week at most in order to get our bachelors degree :)

Mehmet Salih YILDIRIM
-------------- next part --------------
An HTML attachment was scrubbed...
URL: https://software.sandia.gov/pipermail/trilinos-users/attachments/20100527/e60c5f9e/attachment.html