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

$\endgroup$
  • $\begingroup$ Are $A$ and $B$ dense or sparse? $\endgroup$ – Dan May 23 '12 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 May 23 '12 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 May 24 '12 at 3:30
  • $\begingroup$ @Costis precondition here is for iterative eigenvalue solver. $\endgroup$ – bobye May 24 '12 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 May 24 '12 at 13:27
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.

$\endgroup$
  • $\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 May 24 '12 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 May 24 '12 at 14:47
  • $\begingroup$ @bobye: I updated my answer to account for the singular case without any approximation. $\endgroup$ – Arnold Neumaier May 25 '12 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 May 26 '12 at 2:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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