Complex SLEPc gives different eigenvectors

dolfinx
1

Dear All,

I solved very simple eigenvalue in dolfin and dolfinx. I need to use complex build of slepc in dolfinx for future projects. Hence, I am switching over. I have a dolfin code which gives the eigenvectors I want;

import dolfin as dolf
from math import pi
import numpy as np
from slepc4py import SLEPc

mesh = dolf.UnitIntervalMesh(5)
dx = dolf.Measure('dx', domain=mesh)

V = dolf.FunctionSpace(mesh, 'CG', 1)

u = dolf.TrialFunction(V)
v = dolf.TestFunction(V)

a_ = dolf.dot(dolf.grad(u), dolf.grad(v)) * dx
c_ = u * v * dx

A = dolf.PETScMatrix()
dolf.assemble(a_, tensor=A)
A = dolf.as_backend_type(A).mat()

C = dolf.PETScMatrix()
dolf.assemble(c_, tensor=C)
C = dolf.as_backend_type(C).mat()

E = SLEPc.EPS().create()

E.setOperators(A, C)
E.setTarget(np.pi)
E.setWhichEigenpairs(SLEPc.EPS.Which.TARGET_MAGNITUDE)  # TARGET_REAL or TARGET_IMAGINARY
E.setFromOptions()
E.solve()

vr, vi = A.getVecs()
omega = E.getEigenpair(1, vr, vi)

print(omega,vr[:],vi[:])

Output of this code is;

(10.198390006583926+0j) 
[ 0.53452248  0.43243777  0.16517653 -0.16517653 -0.43243777 -0.53452248]
 [0. 0. 0. 0. 0. 0.]

This is expected output since we have real eigenvalue and then zero imaginary part of eigenvector.

And I have the code which does the same thing in dolfinx with complex PETSc and SLEPc;

from dolfinx import Function, FunctionSpace, UnitIntervalMesh
from dolfinx.fem.assemble import assemble_matrix
from mpi4py import MPI
from ufl import Measure, FacetNormal, TestFunction, TrialFunction, dx, grad, inner
import numpy as np
from slepc4py import SLEPc

mesh = UnitIntervalMesh(MPI.COMM_WORLD, 5)

dx = Measure('dx', domain=mesh)

V = FunctionSpace(mesh, ('CG', 1))

u = TrialFunction(V)
v = TestFunction(V)

a  = inner(grad(u), grad(v)) * dx
c  = inner(u,v) * dx

A = assemble_matrix(a)
A.assemble()

C = assemble_matrix(c)
C.assemble()

E = SLEPc.EPS().create()
E.setOperators(A, C)
E.setTarget(np.pi)
E.setWhichEigenpairs(SLEPc.EPS.Which.TARGET_MAGNITUDE)  # TARGET_REAL or TARGET_IMAGINARY
E.setFromOptions()
E.solve()

vr, vi = A.getVecs()
omega = E.getEigenpair(1, vr, vi)

print(omega,vr[:],vi[:])

However, dolfinx script gives eigenvector values which has nonzero complex parts in its entries even if it has same eigenvalue with dolfin script;

(10.19839000658394-2.090963493473814e-15j) 
[ 0.22892832+0.36687089j  0.28297096+0.45347736j  0.08744284+0.14013221j
 -0.08744284-0.14013221j -0.22892832-0.36687089j -0.28297096-0.45347736j]
[0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]

I cannot figure out why they give different eigenvectors? Is there anybody who can say what is going on under the hood for complex SLEPc and PETSc builds?

Thanks in advance.

2

It looks like what’s happening here is that the complex build is choosing a different arbitrary scalar in the eigenvector. Clearly, if \mathbf{A}\mathbf{v} = \lambda\mathbf{v}, then \mathbf{A}(\alpha\mathbf{v}) = \lambda(\alpha\mathbf{v}). It’s just somewhat obscured in this case, because the degrees of freedom appear to be ordered differently in the vector (which I’m guessing is some sort of dolfin vs dolfinx difference in the construction of the FunctionSpace). You can guess the permutation by looking at magnitude, and look at ratios of corresponding entries in the vectors, to see that the dolfinx vector is scaled by a constant with complex modulus one, as illustrated by the following Octave session:

octave:5> 0.53452248 / (0.28297096+0.45347736j)
ans =  0.5294 - 0.8484i
octave:6> 0.43243777 / (0.22892832+0.36687089j)
ans =  0.5294 - 0.8484i
octave:7> 0.16517653 / (0.08744284+0.14013221j)
ans =  0.5294 - 0.8484i
octave:8> 0.5294^2 + 0.8484^2
ans = 1.0000
1 Like
4

Thanks for your answer.

What would happen if I solve polynomial eigenvalue problem? Since have this script;

import sys, slepc4py
slepc4py.init(sys.argv)

from petsc4py import PETSc
from slepc4py import SLEPc

Print = PETSc.Sys.Print


def construct_operators(m,n):
    """
    Standard symmetric eigenproblem corresponding to the
    Laplacian operator in 2 dimensions.
    """
    Print("Quadratic Eigenproblem, N=%d (%dx%d grid)"% (m*n, m, n))
    # K is the 2-D Laplacian
    K = PETSc.Mat().create()
    K.setSizes([n*m, n*m])
    K.setFromOptions( )
    K.setUp()
    Istart, Iend = K.getOwnershipRange()
    for I in range(Istart,Iend):
        v = -1.0; i = I//n; j = I-i*n;
        if i>0:
            J=I-n; K[I,J] = v
        if i<m-1:
            J=I+n; K[I,J] = v
        if j>0:
            J=I-1; K[I,J] = v
        if j<n-1:
            J=I+1; K[I,J] = v
        v=4.0; K[I,I] = v
    K.assemble()
    # C is the zero matrix
    C = PETSc.Mat().create()
    C.setSizes([n*m, n*m])
    C.setFromOptions( )
    C.setUp()
    C.assemble()
    # M is the identity matrix
    M = PETSc.Mat().create()
    M.setSizes([n*m, n*m])
    M.setFromOptions( )
    M.setUp()
    M.assemble()
    M.shift(1)
    #
    return M, C, K

def solve_eigensystem(M, C, K):
    # Setup the eigensolver
    Q = SLEPc.PEP().create()
    A = [ ]
    A.append(K)
    A.append(C)
    A.append(M)
    Q.setOperators(A)
    Q.setDimensions(6)
    Q.setProblemType(SLEPc.PEP.ProblemType.GENERAL)
    Q.setFromOptions()
    # Solve the eigensystem
    Q.solve()
    # Create the result vectors
    xr, xi = K.createVecs()

    its = Q.getIterationNumber()
    Print("Number of iterations of the method: %i" % its)
    sol_type = Q.getType()
    Print("Solution method: %s" % sol_type)
    nev, ncv, mpd = Q.getDimensions()
    Print("")
    Print("Number of requested eigenvalues: %i" % nev)
    tol, maxit = Q.getTolerances()
    Print("Stopping condition: tol=%.4g, maxit=%d" % (tol, maxit))
    nconv = Q.getConverged()
    Print("Number of converged approximate eigenpairs: %d" % nconv)
    if nconv > 0:
        Print("")
        Print("          k           ||(k^2M+Ck+K)x||/||kx|| ")
        Print("-------------------- -------------------------")
        for i in range(nconv):
            k = Q.getEigenpair(i, xr, xi)
            error = Q.computeError(i)
            if k.imag != 0.0: 
                Print("%9f%+9f j    %12g" % (k.real, k.imag, error))
            else: 
                Print("%12f         %12g" % (k.real, error))
    Print("")
    print("XR: ", xr[:])
    print("XI: ", xi[:])

if __name__ == '__main__':
    opts = PETSc.Options()
    m = opts.getInt('m', 4)
    n = opts.getInt('n', m)
    M, C, K = construct_operators(m,n)
    solve_eigensystem(M, C, K)
    M = C = K = None

Dolfin gives me;

XR:  [-0.29835467 -0.28198379 -0.03382559  0.10317374 -0.38515753 -0.29835467
  0.10317374  0.26452908 -0.26452908 -0.10317374  0.29835467  0.38515753
 -0.10317374  0.03382559  0.28198379  0.29835467]
XI:  [ 0.00602335  0.01544338  0.02248575  0.01741814 -0.00197476  0.00602335
  0.01741814  0.0164624  -0.0164624  -0.01741814 -0.00602335  0.00197476
 -0.01741814 -0.02248575 -0.01544338 -0.00602335]

While dolfinx gives;

XR:  [ 0.16139385-0.25541378j  0.22354849-0.30347618j  0.17708191-0.16776581j
  0.08620935-0.03582979j  0.13733913-0.26764639j  0.16139385-0.25541378j
  0.08620935-0.03582979j  0.01568806+0.08764797j -0.01568806-0.08764797j
 -0.08620935+0.03582979j -0.16139385+0.25541378j -0.13733913+0.26764639j
 -0.08620935+0.03582979j -0.17708191+0.16776581j -0.22354849+0.30347618j
 -0.16139385+0.25541378j]
XI:  [0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j
 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]

Then I tried to get;

>>> (-0.29835467 + 0.00602335j)/(0.16139385-0.25541378j)
(-0.5443556463593747-0.8241490199347193j)
>>> (-0.28198379+0.01544338j)/(0.22354849-0.30347618j)
(-0.47668704795634964-0.5780391733769698j)

They are quite different?

5

Since this script doesn’t appear to use dolfin or dolfinx, it’s not clear to me how to get the two sets of output vectors. However, one thing to keep in mind is that the order of degrees of freedom in the vectors is different between dolfin and dolfinx. Notice that, in my previous response, I took the ratio of the 0^\text{th} entry of the dolfin result and the 1^\text{st} entry of the dolfinx result, and the ratio of the the 1^\text{st} entry of the dolfin result and the 0^\text{th} entry of the dolfinx result; I didn’t simply compare entries with the same index. Since there were only a few entries, I just figured out the permutation by guess-and-check; if you want to compare entries in longer vectors, you’ll need to reverse-engineer the mapping between degrees of freedom in dolfin and dolfinx vectors.

1 Like
6

My eigenvectors are quite big, so it is difficult to get mapping. I think the unclearity comes from SLEPc. I will open an issue and will see how it goes.

7

SLEPc team explains this in their FAQ page:

16. Eigenvectors have nonzero imaginary part
A real symmetric matrix has real eigenvectors, but when building SLEPc with complex scalars the computed eigenvectors have nonzero imaginary part. The rationale is the following. In real scalars, if x is a unit-norm eigenvector then -x is also a valid eigenvector. In complex scalars, if x is a unit-norm eigenvector then alpha*x is also a valid eigenvector, where alpha is a generalized sign, i.e., alpha=exp(theta*j) for any theta. So if one wants the imaginary part to be zero, the eigenvectors returned by SLEPc must be normalized a posteriori, as is done for example in ex20.c (or the equivalent python example ex7.py). SLEPc does not know if the input matrix is real or complex, so it cannot normalize the vectors internally.

you find this text here:
https://slepc.upv.es/documentation/faq.htm#:~:text=nonzero%20imaginary%20part-,A%20real%20symmetric,-matrix%20has%20real

Then they propose a solution in ex20.c:

394: PetscErrorCode FixSign(Vec x)
395: {
396:   PetscMPIInt       rank;
397:   PetscScalar       sign=0.0;
398:   const PetscScalar *xx;

400:   PetscFunctionBeginUser;
401:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank));
402:   if (!rank) {
403:     PetscCall(VecGetArrayRead(x,&xx));
404:     sign = *xx/PetscAbsScalar(*xx);
405:     PetscCall(VecRestoreArrayRead(x,&xx));
406:   }
407:   PetscCallMPI(MPI_Bcast(&sign,1,MPIU_SCALAR,0,PETSC_COMM_WORLD));
408:   PetscCall(VecScale(x,1.0/sign));
409:   PetscFunctionReturn(PETSC_SUCCESS);
410: }

I did something like this in python:

#...other code here...

xr, xi      = K.createVecs()

eig_vectors = []

    for i in range(nconv):
        k    = problem.getEigenpair(i, xr, xi)
        vect = xr.getArray()
        sign = vect[0]/np.abs(vect[0])

        vect_mod = vect/sign
        eig_vectors.append(vect_mod.copy())

Hope that this helps!

1 Like