cond2sp — computes an approximation of the 2-norm condition number of a s.p.d. sparse matrix
[K2, lm, vm, lM, vM] = cond2sp(A, C_ptr [, rtol, itermax, verb])
a real symmetric positive definite sparse matrix
a pointer to a Cholesky factorization (got with taucs_chfact)
(optional) relative tolerance (default 1.e-3) (see details in DESCRIPTION)
(optional) maximum number of iterations in the underlying algorithms (default 30)
(optional) boolean, must be %t for displaying the intermediary results, and %f (default) if you don't want.
estimated 2-norm condition number K2 = ||A||_2 ||A^(-1)||_2 = lM/lm
(real positive scalar) minimum eigenvalue
associated eigenvector
(real positive scalar) maximum eigenvalue
associated eigenvector
This quick and dirty function computes (lM,vM) using the iterative
power method and (lm,vm) with the inverse iterative power method, then
K2 = lM/lm. For each method the iterations are stopped until the following
condition is met :
| (l_new - l_old)/l_new | < rtol
but 4 iterations are nevertheless required and also the iterations are stopped
if itermax is reached (and a warning message is issued). As the matrix is symmetric
this is the rayleigh quotient which gives the estimated eigenvalue at each step
(lambda = v'*A*v). You may called this function with named parameter, for
instance if you want to see the intermediary result without setting yourself the
rtol and itermax parameters you may called this function with the syntax :
[K2, lm, vm, lM, vM] = cond2sp(A , C_ptr, verb=%t )
This function is intended to get an approximation of the 2-norm condition number (K2) and
with the methods used, the precision on the obtained eigenvectors (vM and vm) are generally
not very good. If you look for a smaller residual ||Av - l*v||, you may apply some inverse
power iterations from v0 with the matrix :
B = A - l0*speye(A)
For instance, applied 5 such iterations for (lm,vm) is done with :
l0 = lm ; v0 = vm; // or l0 = lM ; v0 = vM; // to polish (lM,vM)
B = A - l0*speye(A);
LUp = umf_lufact(B);
vr = v0; nstep = 5;
for i=1:nstep, vr = umf_lusolve(LUp, vr, "Ax=b", B); vr = vr/norm(vr) ; end
umf_ludel(LUp); // if you don't use anymore this factorization
lr = vr'*A*vr;
norm_r0 = norm(A*v0 - l0*v0);
norm_rr = norm(A*vr - lr*vr);
// Bauer-Fike error bound...
mprintf(" first estimated eigenvalue : l0 = %e \n\t", l0)
mprintf(" |l-l0| <= ||Av0-l0v0|| = %e , |l-l0|/l0 <= %e \n\r", norm_r0, norm_r0/l0)
mprintf(" raffined estimated eigenvalue : lr = %e \n\t", lr)
mprintf(" |l-lr| <= ||Avr-lrvr|| = %e , |l-lr|/lr <= %e \n\r", norm_rr, norm_rr/lr)