[SciPy-Dev] Generalized eigenproblem with rank deficient matrices

[Nils Wagner]

On Sun, 4 Sep 2011 11:12:01 -0600
  Charles R Harris <charlesr.harris at gmail.com> wrote:
> On Sun, Sep 4, 2011 at 10:55 AM, Charles R Harris 
><charlesr.harris at gmail.com
>> wrote:
> 
>>
>>
>> On Sun, Sep 4, 2011 at 10:09 AM, Nils Wagner 
>><nwagner at iam.uni-stuttgart.de
>> > wrote:
>>
>>> On Sun, 4 Sep 2011 09:29:19 -0600
>>>  Charles R Harris <charlesr.harris at gmail.com> wrote:
>>> > On Sun, Sep 4, 2011 at 7:53 AM, Nils Wagner
>>> ><nwagner at iam.uni-stuttgart.de>wrote:
>>> >
>>> >> Hi all,
>>> >>
>>> >> how can I solve the eigenproblem
>>> >>
>>> >> A x = \lambda B x
>>> >>
>>> >> where both matrices are rank deficient ?
>>> >>
>>> >
>>> > I'd do eigh and transform the problem to something 
>>>like:
>>> >
>>> > U * A  * U^t * x= \lambda D * x
>>> >
>>> > where D is diagonal. Note that the solutions may not 
>>>be
>>> >unique and \lambda
>>> > can be arbitrary, as you can see by studying
>>> >
>>> > A = B = array([[1, 0], [0, 0]])
>>> >
>>> > Where there are solutions for arbitrary \lambda.
>>> >Likewise, there may be no
>>> > solutions under the requirement that x is non-zero:
>>> >
>>> > A = array([[1, 1], [1, 0]]),
>>> > B = array([[1, 0], [0, 0]])
>>> >
>>> > The usual case where B is positive definite 
>>>corresponds
>>> >to finding extrema
>>> > on a compact surface x^t * B *x = 1, but the surface 
>>>is
>>> >no longer compact
>>> > when B isn't positive definite. Note that these cases
>>> >are all sensitive to
>>> > roundoff error.
>>> >
>>> > Chuck
>>>
>>> Hi Chuck,
>>>
>>> I am only interested in the real and complex
>>> eigensolutions.
>>> The complex eigenvalues appear in pairs  a_i \pm 
>>>\sqrt{-1}
>>> b_i sind A and B are real.
>>> How can I reject infinite eigenvalues ?
>>> Both matrices, A and B, are indefinite.
>>>
>>>
>> It depends on the particular problem. In general, the 
>>solution to the
>> generalized eigenvalue problem starts by making a 
>>variable substitution that
>> reduces B to the identity matrix, usually by using a 
>>Cholesky factorization,
>> i.e., B = U^t U, y = U x. This can still be done in the 
>>(numerically)
>> indefinite case but Cholesky won't be reliable and that 
>>is why I suggested
>> eigh. Note that the problem we are trying to solve is
>> finding the extrema of x^t A x subject to the constraint 
>>x^t B x = 1,
>> \lambda is then a Lagrange  multiplier. I assume A and B 
>>are both symmetric?
>> Anyway, in terms of y, the problem then reduces to 
>>finding extrema of y^t
>> U^t^{-1} A U^{-1} y subject to y^t D y = 1, where D is 
>>diagonal and has ones
>> along part of the diagonal, zeros for the remainder. In 
>>the usual case, D is
>> the identity. The trick is then to divide y into two 
>>parts, one for where D
>> is one (u), another for the rest (v), so that y = [u v]. 
>>If you are lucky,
>> the v can be solved in terms of u using the transformed 
>>A, and things will
>> reduce to an eigenvalue problem for u. If the original A 
>>was symmetric, so
>> will be the reduced problem. If you can't solve v as a 
>>function of u, then
>> you can reduce things further, but it is possible that 
>>at some point there
>> is no solution.
>>
>> I don't have practical experience with this sort of 
>>problem with indefinite
>> B, so I can't tell you much more than that. I assume 
>>you've googled for
>> relevant documents.
>>
>> It occurs to me that even if A and B aren't symmetric, 
>>that you can still
> bring B to the desired form using the svd.
> 
> Chuck

  
Thank you very much. A and B are symmetric.

Nils