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.

Improve this question
$\endgroup$
3
  • 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, 2013 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, 2013 at 9:23
  • $\begingroup$ Ah, I see now that the $K_i$s are matrices. $\endgroup$
    – Nico Schlömer
    May 14, 2013 at 9:43

1 Answer 1

Reset to default
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.)

Improve this answer
$\endgroup$
26
  • $\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, 2013 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, 2013 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, 2013 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, 2013 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, 2013 at 13:00

Your Answer

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

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