[Numpy-discussion] Getting non-normalized eigenvectors from generalized eigenvalue solution?

[Fahreddın Basegmez]

Howdy,

Is it possible to get non-normalized eigenvectors from scipy.linalg.eig(a, b)?  Preferably just by using  numpy.

BTW, Matlab/Octave provides this with its eig(a, b) function but I would like to use numpy for obvious reasons.

Regards,

Fahri 

[Olivier Delalleau]

I'm probably missing something, but... Why would you want non-normalized eigenvectors?

-=- Olivier

2011/12/20 Fahreddın Basegmez <mang...@gmail.com>

[Fahreddın Basegmez]

I am computing normal-mode frequency response of a mass-spring system.  The algorithm I am using requires it.

_______________________________________________
NumPy-Discussion mailing list
NumPy-Di...@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

[Olivier Delalleau]

Hmm... ok ;) (sorry, I can't follow you there)

Anyway, what kind of non-normalization are you after? I looked at the doc for Matlab and it just says eigenvectors are not normalized, without additional details... so it looks like it could be anything.

[Fahreddın Basegmez]

If I can get the same response as Matlab I would be all set.  

Octave results

>> STIFM

STIFM =

Diagonal Matrix

     1020        0        0        0        0        0

        0     1020        0        0        0        0

        0        0     1020        0        0        0

        0        0        0   102000        0        0

        0        0        0        0   102000        0

        0        0        0        0        0   204000

>> MASSM

MASSM =

Diagonal Matrix

    0.25907          0          0          0          0          0

          0    0.25907          0          0          0          0

          0          0    0.25907          0          0          0

          0          0          0   26.00000          0          0

          0          0          0          0   26.00000          0

          0          0          0          0          0   26.00000

>> [a, b] = eig(STIFM, MASSM)

a =

   0.00000   0.00000   0.00000   1.96468   0.00000   0.00000

   0.00000   0.00000   0.00000   0.00000   1.96468   0.00000

   0.00000   0.00000   1.96468   0.00000   0.00000   0.00000

   0.19612   0.00000   0.00000   0.00000   0.00000   0.00000

   0.00000   0.19612   0.00000   0.00000   0.00000   0.00000

   0.00000   0.00000   0.00000   0.00000   0.00000   0.19612

b =

Diagonal Matrix

   3923.1        0        0        0        0        0

        0   3923.1        0        0        0        0

        0        0   3937.2        0        0        0

        0        0        0   3937.2        0        0

        0        0        0        0   3937.2        0

        0        0        0        0        0   7846.2

Numpy Results

>>> STIFM

array([[   1020.,       0.,       0.,       0.,       0.,       0.],

       [      0.,    1020.,       0.,       0.,       0.,       0.],

       [      0.,       0.,    1020.,       0.,       0.,       0.],

       [      0.,       0.,       0.,  102000.,       0.,       0.],

       [      0.,       0.,       0.,       0.,  102000.,       0.],

       [      0.,       0.,       0.,       0.,       0.,  204000.]])

>>> MASSM

array([[  0.25907,   0.     ,   0.     ,   0.     ,   0.     ,   0.     ],

       [  0.     ,   0.25907,   0.     ,   0.     ,   0.     ,   0.     ],

       [  0.     ,   0.     ,   0.25907,   0.     ,   0.     ,   0.     ],

       [  0.     ,   0.     ,   0.     ,  26.     ,   0.     ,   0.     ],

       [  0.     ,   0.     ,   0.     ,   0.     ,  26.     ,   0.     ],

       [  0.     ,   0.     ,   0.     ,   0.     ,   0.     ,  26.     ]])

>>> a, b = linalg.eig(dot( linalg.pinv(MASSM), STIFM))

>>> a

array([ 3937.15984097,  3937.15984097,  3937.15984097,  3923.07692308,

        3923.07692308,  7846.15384615])

>>> b

array([[ 1.,  0.,  0.,  0.,  0.,  0.],

       [ 0.,  1.,  0.,  0.,  0.,  0.],

       [ 0.,  0.,  1.,  0.,  0.,  0.],

       [ 0.,  0.,  0.,  1.,  0.,  0.],

       [ 0.,  0.,  0.,  0.,  1.,  0.],

       [ 0.,  0.,  0.,  0.,  0.,  1.]])

[Fahreddın Basegmez]

I should include the scipy response too I guess.

scipy.linalg.eig(STIFM, MASSM)

(array([ 3937.15984097+0.j,  3937.15984097+0.j,  3937.15984097+0.j,

        3923.07692308+0.j,  3923.07692308+0.j,  7846.15384615+0.j]), array([[ 1.,  0.,  0.,  0.,  0.,  0.],

       [ 0.,  1.,  0.,  0.,  0.,  0.],

       [ 0.,  0.,  1.,  0.,  0.,  0.],

       [ 0.,  0.,  0.,  1.,  0.,  0.],

       [ 0.,  0.,  0.,  0.,  1.,  0.],

       [ 0.,  0.,  0.,  0.,  0.,  1.]]))

[Olivier Delalleau]

Hmm, sorry, I don't see any obvious logic that would explain how Octave obtains this result, although of course there is probably some logic...

Anyway, since you seem to know what you want, can't you obtain the same result by doing whatever un-normalizing operation you are after?

[Fahreddın Basegmez]

I don't think I can do that.  I can go to the normalized results but not the other way.

[Olivier Delalleau]

What I don't get is that "un-normalized" eigenvectors can be pretty much anything. If you care about the specific output of Matlab / Octave, it means you understand the particular "un-normalization" that these programs use. In that case you should be able to recover it from the normalized output from numpy.

[Fahreddın Basegmez]

I think I am interested in the non-normalized eigenvectors not the un-normalized ones.  Once the eig function computes the generalized eigenvectors I would like to use them as they are. 

I would think this would be a common request since the normal-mode frequency response is used in many different fields like chemical and biomolecular sciences as well as engineering and physics.  Mathematically there may be no difference between the normalized and non-normalized eigenvectors but physically there is.  In my case those values represent deflections.  Advantage of the normal-modes is you can apply damping in each direction independent of each other.  Amount of damping we apply may be dependent on those deflections so I would need to use the non-normalized results.

[Olivier Delalleau]

Ok well I'm sorry, I have no idea what would be the difference between "non-normalized" and "un-normalized".

In PCA you may decide to scale your eigenvectors by the inverse of the square root of their corresponding eigenvalue so that your projected data has unit variance, but it doesn't seem to be what you're after.

Can you point to a link that explains what are the "non-normalized eigenvectors" in your application?

[Fahreddın Basegmez]

Sorry about that.  I don't think that terminology is commonly used.  This is what I mean.

Let's say I solve the equations and compute the eigenvalues and eigenvectors for the given two matrices.  I call these results "non-normalized". Then they can be normalized.  Once they are normalized if I multiply them by any scalar they would become "un-normalized".  They would still be eigenvectors but not necessarily the "non-normalized" ones.  I think there is a specific algorithm that Matlab uses to solve them but I do not know what that algorithm is.

[Lennart Fricke]

Dear Fahreddın,
I think, the norm of the eigenvectors corresponds to some generic
amplitude. But that is something you cannot extract from the solution of
the eigenvalue problem but it depends on the initial deflection or
velocities.

So I think you should be able to use the normalized values just as well
as the non-, un- or not normalized ones.

Octave seems to normalize that way that, transpose(Z).B.Z=I, where Z is
the matrix of eigenvectors, B is matrix B of the generalized eigenvalue
problem and I is the identity. It uses lapack functions. But that's only
true if A,B are symmetric. If not it normalizes the magnitude of largest
element of each eigenvector to 1.

I believe you can get it like that. If U is a Matrix with normalization
factors it is diagonal and Z.A contains the normalized column vectors.
then it is:

transpose(Z.A).B.Z.A
=transpose(A).transpose(Z).B.Z.A
=A.transpose(Z).B.Z.A=I

and thus invert(A).invert(A)=transpose(Z).B.Z
As A is diagonal invert(A) has the reciprocal elements on the diagonal.
So you can easily extract them

A=diag(1/sqrt(diag(transpose(Z).B.Z)))

I hope that's correct.

Best Regards
Lennart

[Olivier Delalleau]

Aaah, thanks a lot Lennart, I knew there had to be some logic to Octave's output, but I couldn't see it...

-=- Olivier

2011/12/21 Lennart Fricke <pge0...@studserv.uni-leipzig.de>

[Andrew Jaffe]

Just to be completely clear, there is no such thing as a
"non-normalized" eigenvector. An eigenvector is only determined *up to a
scalar normalization*, which is obvious from the eigenvalue equation:

A v = l v

where A is the matrix, l is the eigenvalue, and v is the eigenvector.
Obviously v is only determined up to a constant factor. A given eigen
routine can return anything at all, but there is no native
"non-normalized" version.

Traditionally, you can decide to return "normalized" eigenvectors with
the scalar factor determined by norm(v)=1 for some suitable norm. (I
could imagine that an algorithm could depend on that.)

Andrew

On 21/12/2011 07:01, Olivier Delalleau wrote:
> Aaah, thanks a lot Lennart, I knew there had to be some logic to
> Octave's output, but I couldn't see it...
>
> -=- Olivier
>
> 2011/12/21 Lennart Fricke <pge0...@studserv.uni-leipzig.de

> <mailto:pge0...@studserv.uni-leipzig.de>>

> NumPy-Di...@scipy.org <mailto:NumPy-Di...@scipy.org>
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

[Fahreddın Basegmez]

According to this page eigenvectors are normalized with respect to the second matrix.  Do you guys have any idea how that's done?

http://www.kxcad.net/Altair/HyperWorks/oshelp/frequency_response_analysis.htm

 "If the eigenvectors are normalized with respect to the mass matrix, the modal mass matrix is the unity matrix and the modal stiffness matrix is a diagonal matrix holding the eigenvalues of the system. This way, the system equation is reduced to a set of uncoupled equations for the components of d that can be solved easily."

[Lennart Fricke]

I read it like that:
(**T is the transpose)

Let's call M the mass matrix and N the modal mass matrix. Then
X**T*M*X=N. If X (matrix of eigenvectors) is normalized with respect to
M, N is I (unity) so it just mean that X**T*M*X=I. That is what octave
and matlab give you.

For this to be true. x**T*M*x=1 must be true for each column x of X.
Thus if y is the not normalized eigenvector.
a*y**T*M*a*y=1
a**2 * y**T*M*y=1
a**2 = 1/(y**T*M*y)

For me the question, why X diagonalizes M and K at the same time, remains.

Another hint. If the matrizes are hermitian (symmetric if real) you can
use scipy.linalg.eigh, which gives you result with correct
normalization.

Best regards
Lennart

_______________________________________________
NumPy-Discussion mailing list
NumPy-Di...@scipy.org

http://mail.scipy.org/mailman/listinfo/numpy-discussion