How to compute the condition number efficiently?

Linear Algebra
1

How to compute the condition number of a problem efficiently? I know one can convert the matrix to dense and use numpy, but I am looking for something more efficient.

from fenics import *
import numpy as np

mesh = UnitIntervalMesh(10)
V = FunctionSpace(mesh, "Lagrange", 1)
u = TrialFunction(V)
v = TestFunction(V)
F = inner(v, u) * dx
A = assemble(lhs(F))
np.linalg.cond(A.array())
2

Consider using the SLEPcEigenSolver:

from dolfin import *
import numpy as np
import time

mesh = UnitSquareMesh(60,60)
V = FunctionSpace(mesh, "CG", 1)
u,v = TrialFunction(V), TestFunction(V)
a = inner(u,v)*dx

A = PETScMatrix()
assemble(a, tensor=A)


start = time.time()
eigenSolver = SLEPcEigenSolver(A)
eigenSolver.parameters["spectrum"]="smallest magnitude"
eigenSolver.solve(1)
eigen_min = eigenSolver.get_eigenvalue(0)[0]

eigenSolver.parameters["spectrum"]="largest magnitude"
eigenSolver.solve(1)
eigen_max = eigenSolver.get_eigenvalue(0)[0]
end = time.time()
print("Condition number {0:.2e} ({1:.2e} s)".format(eigen_max/eigen_min,
                                                    end-start))

start = time.time()
cond_numpy = np.linalg.cond(A.array())
end = time.time()
print("Numpy condition number {0:.2e} ({1:.2e} s)".format(cond_numpy,
                                                          end-start))

Note that you might have to increase the argument in

eigenSolver.solve(1)

when increasing the number of cells in your mesh.

2 Likes
3

Nice! I was not patient enough to wait for the numpy result in your example :smiley: . Is there some docs on the eigen solver? E.g. the following questions pop up:

  • What does the 1 in eigenSolver.solve(1) mean?
  • What is eigenSolver.parameters["spectrum"] about?
  • What do the zeros in eigenSolver.get_eigenvalue(0)[0] mean?

Also where can I read about PETScMatrix? e.g. what sparsity structure does it have and how do I linear algebra with it?

4

See for instance this demo for information about SLEPcEigen Solvers.

PETScMatrix is just a wrapper around PETSc’s sparse matrices and its docs can be found here.

1 Like
5

Also how to add a boundary condition before computing cond?

from fenics import *
import numpy as np

mesh = UnitIntervalMesh(10)
V = FunctionSpace(mesh, "Lagrange", 1)
u = TrialFunction(V)
v = TestFunction(V)
F = inner(v, u) * dx
A = assemble(lhs(F))
bc = DirichletBC(V, 1, "on_boundary")
bc.apply(A)
np.linalg.cond(A.array())
6

Please see the demo I linked you to.

1 Like
7

The SLEPcEigenSolver computes a partial spectrum. The “spetrum” argument lets you choose, which part of the spectrum should be computed. The eigenSolver.solve(1) means that roughly one eigenvalue should be computed. re, im = solver.get_eigenvalue(i) gives the real and imaginary part of the ith eigenvalue.

8 

def compute_eig0(solver, neig_min=1, neig_max=1000):
    neig = neig_min
    ret = None
    while True:
        neig = min(neig, neig_max)
        solver.solve(neig)
        try:
            ret = solver.get_eigenvalue(0)
            break
        except Exception:
            pass
        if neig >= neig_max:
            raise f"Reached neig_max = {neig_max}"
        neig = 2*neig
    return ret

def cond(F, bcs):
    A = PETScMatrix()
    dummy = rhs(F)
    assemble_system(lhs(F), dummy, A_tensor=A, bcs=bcs)
    eigenSolver = SLEPcEigenSolver(A)
    eigenSolver.parameters["spectrum"]="smallest magnitude"
    eigen_min, _ = compute_eig0(eigenSolver)
    
    eigenSolver.parameters["spectrum"]="largest magnitude"
    eigen_max, _ = compute_eig0(eigenSolver)
    return abs(eigen_max / eigen_min)