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