Compute Null Space of Assembled PETSc Matrix?

Linear Algebra
1

I have rank-deficient PETSc Matrix that is assembled from a variational problem. I would like to extract the basis of the null space of this matrix. A dense SVD computation works well up to a certain problem size. However, the problems I’d eventually like to solve are much larger and would benefit from parallelization.

Some approaches I’ve tried:

  • numpy SVD: a dense solver that works quickly until matrices in the ~4000x4000 range (what I’m currently using)

  • Scipy null_space(): operates in a dense mode and doesn’t seem to scale well.

  • Scipy sparse SVD, svds(): Convergence for smallest singular vectors is extremely slow or never converges with arpack or lobpcg solvers.

  • SLEPc SVD: Multiple solvers with mixed results. Code snippets below.

    • The ‘debugging’ lapack solver is a dense solver just like numpy and also doesn’t run in parallel:
#SOME VARIATIONAL FORM ASSEMBLED PREVIOUSLY
A = assemble_matrix(form(FORM))
A.assemble()
import slepc4py,sys
slepc4py.init(sys.argv)
Prob = slepc4py.SLEPc.SVD()
Prob.create()
Prob.setOperators(A)
Prob.setType("lapack")
Prob.setTolerances(1E-8, 100)
Prob.solve()
  • The trlanczos solver. I’ve been able to converge some results to reasonable tolerances. but they require massive numbers of iterations and perform significantly worse than the dense solvers for small problems:
#SOME VARIATIONAL FORM ASSEMBLED PREVIOUSLY
A = assemble_matrix(form(FORM))
A.assemble()
Prob = slepc4py.SLEPc.SVD().create()
Prob.setOperators(A)
Prob.setType("trlanczos")
Prob.setDimensions(12,PETSc.DECIDE,PETSc.DECIDE)
Prob.setWhichSingularTriplets(Prob.Which.SMALLEST)
Prob.setTolerances(1e-8, 10000)
Prob.setImplicitTranspose(True)
Prob.solve()

Since the conditioning of the assembled matrix is extremely poor (effectively infinite), convergence is slow for the latter singular vectors (the ones I want to extract). This problem seems to get worse and worse as the problem gets larger and larger. I can provide more information on the structure of my problem if need be, but I’m curious if there is more general advice on how to approach a problem

Is there any best practice to extract the nullspace basis from a rank-deficient assembled PETSc matrix? Or are there any other approaches that I’m not considering above?

1 Like
2

Hi jkrokowski,

What type of variational problem are you solving? If you know that the system matrix you assemble based on the variational form is rank-deficient, then one could determine the nullspace of the system matrix from looking at what types of functions that trivially solve the variational problem.

For example: for the Poisson problem

\qquad \qquad \qquad -\nabla^2u = f \quad \mathrm{in} \ \Omega

with pure Neumann boundary conditions

\qquad \qquad \qquad \nabla u\cdot\mathbf{n} = 0 \quad \mathrm{on} \ \partial\Omega

on the whole boundary \partial\Omega of the domain \Omega, any constant function would satisfy the variational form. From this, we can deduce that a vector of constants will be the nullspace vector of the assembled system matrix.

Let me know if I interpreted your question correctly, and if this is of help.

Cheers,
Halvor

2 Likes