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?