6
$\begingroup$

I need to solve a real generalized eigenvalue problem

$Ax= \lambda Bx(*)$

A and B are calculated from equations below:

$$A=\sum_{i,j=1}^{N}W_{ij}(K_{i}-K_{j})\beta\beta^{T}(K_{i}-K_{j})^{T}$$

$$B=\sum_{i=1}^{N}D_{ii}K_{i}\beta\beta^{T}K_{i}^{T}$$.

where $W$ is a real symmetric $N*N$ matrix with diagonal entries being $0$ and off-diagonal entries between $(0,1)$.

$D$ is an $N*N$ diagonal matrix with $D_{ii}=\sum_{j=1}^NW_{ij}$.

$K_i$ is an $N*M$ matrix with all entries positive.

$\beta>0$ is an $M$ dimensional column vector.

From above equations, A and B should be symmetric semi-definite and B should be positive definite(I did some proof myself).

Maybe because some numerical losses( I are not sure :( ), $B$ appears to have small negative eigenvalues( I do the eigenvalue decomposition using LAPACK routine dsyev() ) and $(*)$ gives complex eigenvalues.

I want to select P smallest eigenvalues of this generalized eigenvalue problem, so complex values here are really a problem. Is there any way to avoid complex eigenvalues in such a case?

By the way I used armadillo as linear algebra library and solve $(*)$ directly using LAPACK routine dggev().

Any suggestions will be appreciated.

$\endgroup$
  • 1
    $\begingroup$ I assume there's a mistake in how $A$ and $B$ are defined: Right now, they both appear to be rank-1-matrices, $A=\alpha_1 \beta \beta^T$, $B=\alpha_2 \beta \beta^T$. $\endgroup$ – Nico Schlömer May 14 '13 at 9:12
  • $\begingroup$ Hmm I guess the rank of the sum of $N$ rank-1 matrices is not always rank-1? $A$and$B$ are $N$ sums of what you said rank-1 matrices. $\endgroup$ – ZeyuHu May 14 '13 at 9:23
  • $\begingroup$ Ah, I see now that the $K_i$s are matrices. $\endgroup$ – Nico Schlömer May 14 '13 at 9:43
5
$\begingroup$

If, as you say, you are sure that you have a symmetric-definite pencil (that is, $\mathbf A$ is symmetric, and $\mathbf B$ is symmetric positive-definite), then LAPACK already has something for directly handling your problem: dsygv(). What it does is to perform a Cholesky decomposition of $\mathbf B$ (if in fact your $\mathbf B$ is not symmetric positive-definite, then you should see a warning), after which the Cholesky triangle thus produced is used to convert your generalized eigenproblem into a regular symmetric eigenproblem that can be solved with all the usual methods. Since this method requires the inversion of the Cholesky triangle of $\mathbf B$, you'll probably want to check if $\mathbf B$ is well-conditioned; you can use dsycon() for the purpose.


There is an alternative method based on the eigendecomposition of $\mathbf B$ if the Cholesky route fails, also discussed in Golub and Van Loan. Briefly, the procedure proceeds like so: give the eigendecomposition $\mathbf B=\mathbf V\mathbf D\mathbf V^\top$, form the matrix $\mathbf W=\mathbf V\sqrt{\mathbf D}$, where $\sqrt{\mathbf D}$ is done by taking the square roots of the diagonal elements. Having formed $\mathbf W$, form $\mathbf C=\mathbf W^{-1}\mathbf A\mathbf W^{-\top}$, which has the same eigenvalues as the pencil $(\mathbf A,\mathbf B)$. (I'll leave the procedure of how to form the eigenvectors as an exercise.)

$\endgroup$
  • $\begingroup$ I mean $B$ should be positive definite, because when given an arbitrary $N$ dimensional column vector $x\neq 0$ and let $v_{i}=K_{i}\beta$ which is also an $N$ dimensional column vector and $v_{i}>0$, we have $$x^TBx=\sum_{i=1}^{N}x^Tv_{i}v_{i}^Tx=\sum_{i=1}^{N}(x^Tv_{i})^2>0$$.But numerically it isn't for I get negative eigenvalues of $B$ after dsyev(). I want to find a way to make this $B$ "numerically positive definite" before going any further using what your quoted dsygv(). Could this be possible or when $B$ is bad-conditioned this problem could not be efficiently solved? $\endgroup$ – ZeyuHu May 14 '13 at 12:10
  • $\begingroup$ That's why I told you to check with dsycon(); what is it returning for your right matrix? $\endgroup$ – J. M. May 14 '13 at 12:21
  • $\begingroup$ Ah, I get an error using dsycon(), but I do get reciprocal of the condition number using $RCOND = 1 / (norm() * norm(inv()))$.Here $RCOND(B)=1.2296e-20$, Hmm It's really ill-conditioned right? $\endgroup$ – ZeyuHu May 15 '13 at 2:47
  • $\begingroup$ Then yes, your $\mathbf B$ is certainly ill-conditioned. I'll edit my post later to include an alternative method. Unfortunately, that method isn't explicitly implemented in LAPACK, so you'll have to implement it yourself. $\endgroup$ – J. M. May 15 '13 at 2:49
  • 1
    $\begingroup$ I think I know how to form the eigenvectors:$Ax=\lambda Bx\Leftrightarrow B^{-1}Ax=\lambda x$,after we do the eigendecompositon $B=VDV^{T}$ and form $B=WW^{T}$, we get an invertible $W$, so $B^{-1}=W^{-T}W^{-1}$, then we have$W^{-T}W^{-1}Ax=\lambda x$, multiply $W^T$ on the left we get $W^{-1}Ax=\lambda W^{T}x$, by adding $W^T$ we form $$W^{-1}AW^{-T}(W^{T}x)=\lambda (W^{T}x)$$, let $C=W^{-1}AW^{-T}$, so we see $C$ has eigenvectors $W^{T}x$, is that correct? But another question came up when doing the eigendecompositon of $B$, I get some small negative doubles, should I simply abandon them? $\endgroup$ – ZeyuHu May 15 '13 at 13:00

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.