[Trilinos-Users] condition number of an Epetra_FECrsMatrix

[Pavel Jiránek]

Hi Nate,

I believe that this is due to the fact that both the Matlab functions
"condest" and "rcond" estimate the condition number of a given matrix with
respect to the 1-norm while the Aztec's estimators (using an eigen- or
singular-value analysis of the coefficient matrix from the Lanczos or
Arnoldi process) estimate (most likely) the 2-norm. Condition numbers with
respect to these norms are equivalent in the sense that one is a scalar
multiple of the other with the scalar lying somewhere in the interval [1/n,
n], where n is the matrix dimension (for square matrices).

In order to verify your results from Aztec, you could try for a small
matrix to use instead the function "cond" on the matrix converted to the
full matrix (cond(full(A))) which computes the 2-norm condition number
using the singular value decomposition. Also note that the quality of the
estimate obtained from the coefficient matrix in Lanczos and Arnoldi
process (used by Aztec) is strongly dependent on the "quality" of the
projection (which depends on the convergence of the extreme eigen/singular
values). You could also use probably some more "advanced" methods for
eigenvalue problems in Anasazi (probably LOBPCG?) if you need an accurate
result.

Best wishes,

Pavel



2013/4/17 Nate Roberts <nate at nateroberts.com>

> Hello all,
>
> I'm interested in studying the condition number of our global system
> matrix, which we store as an Epetra_FECrsMatrix.   My approach thus far
> doesn't give accurate estimates, and I'd like advice on getting better ones.
>
> Based on some sample code I found in the AztecOO documentation, I wrote
> the following method:
>
> double Solution::conditionNumberEstimate( Epetra_LinearProblem & problem
> ) {
>   AztecOOConditionNumber conditionEstimator;
>   conditionEstimator.initialize(*problem.GetOperator());
>
>
>   int maxIters = 40000;
>   double tol = 1e-10;
>   int status = conditionEstimator.computeConditionNumber(maxIters, tol);
>   if (status!=0)
>     cout << "status result from computeConditionNumber(): " << status <<
> endl;
>   double condest = conditionEstimator.getConditionNumber();
>
>
>   return condest;
> }
>
> (The input Epetra_LinearProblem is initialized with the global system
> matrix.)
>
> For some smaller problems, I also tried exporting the matrix to file, and
> loading it into MATLAB, and using rcond() and/or condest() there.  What I
> found was that the method above significantly underreports the condition
> number.  For example, for a fairly dense 243x243 matrix (corresponding to a
> 2x2 finite element mesh), the method above reports 2611.05 as the condition
> number.  MATLAB's condest and rcond both report 16716.  (Decreasing the
> tolerance to 1e-12 and increasing maxIters to 100000, I instead get 2267.47
> from the method above.)
>
> (It may be worth noting that the Epetra_FECrsMatrix ends up storing quite
> a few zeros.  Of the 19237 entries stored, 5024 of them are zeros.)
>
> I'd really like to be able to get high-fidelity condition numbers for this
> matrix without going to MATLAB.  I don't mind if it takes some
> computational effort—I'm not using the condition number during my real
> solves, but studying the effect of various parameters (e.g. choice of basis
> functions) on conditioning.
>
> Any advice would be greatly appreciated!
>
> Thanks,
> Nate
>
> _______________________________________________
> Trilinos-Users mailing list
> Trilinos-Users at software.sandia.gov
> http://software.sandia.gov/mailman/listinfo/trilinos-users
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: https://software.sandia.gov/pipermail/trilinos-users/attachments/20130417/2aa5aff5/attachment.html