[petsc-users] Singular jocabian using SNES

[Marcel Huysegoms]

Hello everyone,

I have a large system of nonlinear equations for which I'm trying to find the optimal solution.
In order to get familiar with the SNES framework, I created a standalone python script (see below), which creates a set of 2D points and transforms them using an affine transformation. The optimizer should then "move" the points back to their original position given the jacobian and the residual vector.

Now I have 2 questions regarding the usage of SNES.

- As in my real application the jacobian often gets singular (contains many rows of only zeros), especially when it approaches the solution. This I simulate below by setting the 10th row equal to zero in the fill-functions. I read (https://scicomp.stackexchange.com/questions/21781/newtons method-goes-to-zero-determinant-jacobian) that quasi-newton approaches like BFGS might be able to deal with such a singular jacobian, however I cannot figure out a combination of solvers that converges in that case.

I always get the message: Nonlinear solve did not converge due to DIVERGED_INNER iterations 0. What can I do in order to make the solver converge (to the least square minimum length solution)? Is there a solver that can deal with such a situation? What do I need to change in the example script?

- In my real application I actually have an overdetermined MxN system. I've read in the manual that the SNES package expects a square jacobian. Is it possible to solve a system having more equations than unknowns?

Many thanks in advance,
Marcel

-----------------------------------------


import sys
import petsc4py
import numpy as np

petsc4py.init(sys.argv)
from petsc4py import PETSc

def fill_function(snes, x, f, points_original):
    x_values = x.getArray(readonly=True)
    diff_vectors = points_original.ravel() - x_values
    f_values = np.square(diff_vectors)
    # f_values[10] = 0
    f.setValues(np.arange(f_values.size), f_values)
    f.assemble()

def fill_jacobian(snes, x, J, P, points_original):
    x_values = x.getArray(readonly=True)
    points_original_flat = points_original.ravel()
    deriv_values = -2*(points_original_flat - x_values)
    # deriv_values[10] = 0
    for i in range(x_values.size):
        P.setValue(i, i, deriv_values[i])
    # print(deriv_values)
    P.assemble()

# ---------------------------------------------------------------------------------------------

if __name__ == '__main__':
    # Initialize original grid points
    grid_dim = 10
    grid_spacing = 100
    num_points = grid_dim * grid_dim
    points_original = np.zeros(shape=(num_points, 2), dtype=np.float64)
    for i in range(grid_dim):
        for j in range(grid_dim):
            points_original[i*grid_dim+j] = (i*grid_spacing, j*grid_spacing)

    # Compute transformed grid points
    affine_mat = np.array([[-0.5, -0.86, 100], [0.86, -0.5, 100]]) # createAffineMatrix(120, 1, 1, 100, 100)
    points_transformed = np.matmul(affine_mat[:2,:2], points_original.T).T + affine_mat[:2,2]

    # Initialize PETSc objects
    num_unknown = points_transformed.size
    mat_shape = (num_unknown, num_unknown)
    A = PETSc.Mat()
    A.createAIJ(size=mat_shape, comm=PETSc.COMM_WORLD)
    A.setUp()
    x, f = A.createVecs()

    options = PETSc.Options()
    options.setValue("-snes_qn_type", "lbfgs")  # broyden/lbfgs
    options.setValue("-snes_qn_scale_type", "none")     # none, diagonal, scalar, jacobian,
    options.setValue("-snes_monitor", "")
    # options.setValue("-snes_view", "")
    options.setValue("-snes_converged_reason", "")
    options.setFromOptions()

    snes = PETSc.SNES()
    snes.create(PETSc.COMM_WORLD)
    snes.setType("qn")
    snes.setFunction(fill_function, f, args=(points_original,))
    snes.setJacobian(fill_jacobian, A, None, args=(points_original,))
    snes_pc = snes.getNPC()     # Inner snes instance (newtonls by default!)
    # snes_pc.setType("ngmres")
    snes.setFromOptions()

    ksp = snes_pc.getKSP()
    ksp.setType("cg")
    ksp.setTolerances(rtol=1e-10, max_it=40000)
    pc = ksp.getPC()
    pc.setType("asm")
    ksp.setFromOptions()

    x.setArray(points_transformed.ravel())
    snes.solve(None, x)


------------------------------------------------------------------------------------------------
------------------------------------------------------------------------------------------------
Forschungszentrum Juelich GmbH
52425 Juelich
Sitz der Gesellschaft: Juelich
Eingetragen im Handelsregister des Amtsgerichts Dueren Nr. HR B 3498
Vorsitzender des Aufsichtsrats: MinDir Volker Rieke
Geschaeftsfuehrung: Prof. Dr.-Ing. Wolfgang Marquardt (Vorsitzender),
Karsten Beneke (stellv. Vorsitzender), Prof. Dr.-Ing. Harald Bolt
------------------------------------------------------------------------------------------------
------------------------------------------------------------------------------------------------

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20210323/6c5c9885/attachment.html>