[petsc-users] Singular jocabian using SNES

[Matthew Knepley]

On Tue, Mar 23, 2021 at 12:39 PM Marcel Huysegoms <m.huysegoms at fz-juelich.de>
wrote:

> 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?
>

SNES is only for solving systems of nonlinear equations. If you want
optimization (least-square, etc.) then you want to formulate your
problem in the TAO interface. It has quasi-Newton methods for those
problems, and other methods as well. That is where I would start.

  Thanks,

     Matt


> Many thanks in advance,
> Marcel
>
> -----------------------------------------
>
> import sysimport petsc4pyimport 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
>
> ------------------------------------------------------------------------------------------------
>
> ------------------------------------------------------------------------------------------------
>
>

-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20210323/fbb89e40/attachment-0001.html>