# generalized eigenvalue problem

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.

• 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$. – Nico Schlömer May 14 '13 at 9:12
• 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. – ZeyuHu May 14 '13 at 9:23
• Ah, I see now that the $K_i$s are matrices. – Nico Schlömer May 14 '13 at 9:43

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.)
• 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? – ZeyuHu May 14 '13 at 12:10
• That's why I told you to check with dsycon(); what is it returning for your right matrix? – J. M. May 14 '13 at 12:21
• 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? – ZeyuHu May 15 '13 at 2:47
• 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. – J. M. May 15 '13 at 2:49
• 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? – ZeyuHu May 15 '13 at 13:00