2

Context: I have some theoretical a posteriori error bounds for a two-grid finite-element scheme that I am using for the solution of the buckling eigenvalue problem. However, one of the terms is prohibitively costly to compute and I wonder if I am doing it in a naive way.

Problem: Letting A denote the my (sparse) stiffness matrix I compute the Cholesky factorisation:

L = chol(A,'lower');

I then need to compute the 2-norm of inv(L). I am currently using 'inv' and svds:

Linv = inv(L);
LinvNorm = svds(Linv,1);

Note that I am using inv as it is in fact the recommended syntax on the Mathworks website: "For sparse inputs, inv(X) creates a sparse identity matrix and uses backslash, X\speye(size(X))."

Question: This is very slow, of course (especially the computation of the inverse, although svds isn't cheap either). Am I missing a trick here?

Edit: I have tried using svds(L,1,0) (to take the reciprocal), but this fails to converge. I am on 2015a so cannot increase the Krylov dimension, which is the recommended fix in newer releases.

Kino
  • 206
  • 2
  • 6
  • 1
    Nothing missing. You are doing an expensive maths operation, it takes time – Ander Biguri Aug 22 '17 at 15:10
  • @AnderBiguri I had thought as much, I was just being optimistic - a better place to ask may have been Mathematics SE. I will have to try to bound these terms, in that case. – Kino Aug 22 '17 at 15:20
  • 1
    @Kino You may also want to ask on Computational Science SE: https://scicomp.stackexchange.com/ – rayryeng Aug 22 '17 at 16:41

1 Answers1

1

What about 

LinvNorm = 1 / svds(L, 1, 0)

i.e. the reciprocal of the smallest singular value of L? (The singular values of the inverse matrix are the reciprocals of the singular values of the original matrix).

Alternative expression are

LinvNorm = 1 / sqrt(eigs(A, 1, 0))

or

LinvNorm = 1 / sqrt(svds(A, 1, 0))

which derive from A == L*L'.

Stefano M
  • 4,267
  • 2
  • 27
  • 45
  • Sorry, yes I did try this initially but svds fails to converge in this case. I am on 2015a so cannot increase the dimension of the Krylov space, which is the recommended fix for this in newer releases. – Kino Aug 23 '17 at 12:24
  • @kino Did you try to compute `eigs(A,1,0)`? Maybe this is more stable, (and for sure more sparse). – Stefano M Aug 23 '17 at 20:18
  • @kino If the above suggestion, based on elementary linear algebra, is not sufficient, you can try to ask a new question on https://scicomp.stackexchange.com/ about a better way of estimating your error bounds. – Stefano M Aug 23 '17 at 20:22