5
$\begingroup$

suppose solving sequential generalized eigenvalue problems

$$A_i x= \lambda Bx, i=1,2,3,\ldots $$

In general setting, we always need to perform LU for matrix B (preconditioned) before to apply the rest iterative algorithm. Is there a numerical library(I have programming experiences with PETSc+SLEPc) or a toolkit that can allow me to separate those two parts, thus to perform LU only once?

By default, LU factorization of $B$ is by direct solver, whose costs may be somewhat comparable, I suppose.

Update: thanks to Arnold, but I want to modify my problem a little, where $B$ has a null vector s.t. $B\mathbf{1}=\mathbf{0},\quad \mathrm{rank}(B) = n-1$ where $A_i,B$ are both $n\times n$ sparse symmetric matrix

Improve this question
$\endgroup$
6
  • $\begingroup$ Are $A$ and $B$ dense or sparse? $\endgroup$
    – Dan
    Commented May 23, 2012 at 22:30
  • $\begingroup$ Are you sure you mean preconditioned? I see no point to preconditioning B if you are using a direct solver to obtain an LU decomposition for it. $\endgroup$
    – Costis
    Commented May 23, 2012 at 23:07
  • 1
    $\begingroup$ Which eigenvalues are you looking for and with which method? If you are looking for small eigenvalues using shift-and-invert, then won't you need a preconditioner for $K_i = A_i - \alpha B$? If that is the case, then no, there is not a reliable way to reuse a factorization of $K_i$ when solving with $K_j$ (though if they are very small or if $K_j$ is well approximated by $K_i + (\text{low rank})$, then there are opportunities. $\endgroup$
    – Jed Brown
    Commented May 24, 2012 at 3:30
  • $\begingroup$ @Costis precondition here is for iterative eigenvalue solver. $\endgroup$
    – bobye
    Commented May 24, 2012 at 13:24
  • $\begingroup$ @JedBrown yeah, I have already invert the problem, and expected to find large eigenvalues. Arnold give a solution when B is nonsingular, but unfortunately B has a null space, $B \mathbf{1} = \mathbf{0}$. I hope to have some way to work around shift transform. $\endgroup$
    – bobye
    Commented May 24, 2012 at 13:27

1 Answer 1

Reset to default
5
$\begingroup$

Factor $PB=LU$ yourself and write a routine for evaluating $L^{-1}PAU^{-1}x$ given $x$ (using two backsolves). Then you can solve the problem $(L^{-1}PAU^{-1}-\lambda I)z=0$ with a standard iterative solver for the ordinary eigenvalue problem.

If $B$ is singular, compute a left null space basis consisting of the rows of $M$, and a right null space basis consisting of the columns of $N$, so that $MB=0$ and $BN=0$. Then you can replace the eigenvalue problem by the modified problem with the matrices $A'=\pmatrix{A & sAN\\MA & sC}$ and $B'=\pmatrix{B & AN\\MA & C}$, with $C$ and $s$ arbitrary. If $x$ solves the original eigenvalue problem then $x'=\pmatrix{x\\0}$ is an eigenvector of the new problem with the same eigenvalue. Now the kernel of the matrix $\pmatrix{B\\MA}$ is trivial since otherwise the eigenvalue problem is ill-posed. This implies that $B'$ is a nonsingular matrix. Thus one can apply the preceding to the modified problem.

The new eigenvalue problem also has the eigenvalue $s$, attained for all vectors of the form $x'=\pmatrix{0\\z}$. Therefore one should choose $s$ such that it lies somewhere in the middle of the expected spectrum.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Very useful, but I found in my situation, $B$ has null space given by $\mathbf{1}$(in fact, stiff matrix derived from FEM), if I apply shift rule $B-\alpha A_i$, I would compute LU each time. Is there a way to work around? $\endgroup$
    – bobye
    Commented May 24, 2012 at 13:21
  • $\begingroup$ You can probably change the near zero entry in the $U$ factor so something of order $\sqrt{\epsilon}$ times the norm of $U$ without seriously degrading performance. $\endgroup$
    – Arnold Neumaier
    Commented May 24, 2012 at 14:47
  • $\begingroup$ @bobye: I updated my answer to account for the singular case without any approximation. $\endgroup$
    – Arnold Neumaier
    Commented May 25, 2012 at 12:35
  • $\begingroup$ many thanks, I think it is a very reasonable approach, that does work for me. I will try to implement. $\endgroup$
    – bobye
    Commented May 26, 2012 at 2:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.