<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
</head>
<body>
Hello everyone,<br>
<br>
I have a large system of nonlinear equations for which I'm trying to find the optimal solution.<br>
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.<br>
<br>
Now I have 2 questions regarding the usage of SNES.<br>
<br>
- 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 (<a class="moz-txt-link-freetext" href="https://scicomp.stackexchange.com/questions/21781/newtons">https://scicomp.stackexchange.com/questions/21781/newtons</a>
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.<br>
<br>
I always get the message: <i>Nonlinear solve did not converge due to DIVERGED_INNER iterations 0.</i> 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?<br>
<br>
- 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?<br>
<br>
Many thanks in advance,<br>
Marcel<br>
<br>
-----------------------------------------<br>
<br>
<pre style="background-color:#f1f1ea;color:#000000;font-family:'DejaVu Sans Mono';font-size:10,5pt;"><span style="color:#000080;font-weight:bold;">import </span>sys
<span style="color:#000080;font-weight:bold;">import </span>petsc4py
<span style="color:#000080;font-weight:bold;">import </span>numpy <span style="color:#000080;font-weight:bold;">as </span>np
petsc4py.init(sys.argv)
<span style="color:#000080;font-weight:bold;">from </span>petsc4py <span style="color:#000080;font-weight:bold;">import </span>PETSc
<span style="color:#000080;font-weight:bold;">def </span>fill_function(<span style="color:#808080;">snes</span>, x, f, points_original):
x_values = x.getArray(<span style="color:#660099;">readonly</span>=<span style="color:#000080;">True</span>)
diff_vectors = points_original.ravel() - x_values
f_values = np.square(diff_vectors)
<span style="color:#808080;font-style:italic;"># f_values[10] = 0
</span><span style="color:#808080;font-style:italic;"> </span>f.setValues(np.arange(f_values.size), f_values)
f.assemble()
<span style="color:#000080;font-weight:bold;">def </span>fill_jacobian(<span style="color:#808080;">snes</span>, x, <span style="color:#808080;">J</span>, P, points_original):
x_values = x.getArray(<span style="color:#660099;">readonly</span>=<span style="color:#000080;">True</span>)
points_original_flat = points_original.ravel()
deriv_values = -<span style="color:#0000ff;">2</span>*(points_original_flat - x_values)
<span style="color:#808080;font-style:italic;"># deriv_values[10] = 0
</span><span style="color:#808080;font-style:italic;"> </span><span style="color:#000080;font-weight:bold;">for </span>i <span style="color:#000080;font-weight:bold;">in </span>range(x_values.size):
P.setValue(i, i, deriv_values[i])
<span style="color:#808080;font-style:italic;"># print(deriv_values)
</span><span style="color:#808080;font-style:italic;"> </span>P.assemble()
<span style="color:#808080;font-style:italic;"># ---------------------------------------------------------------------------------------------
</span><span style="color:#808080;font-style:italic;">
</span><span style="color:#000080;font-weight:bold;">if </span>__name__ == <span style="color:#008000;font-weight:bold;">'__main__'</span>:
<span style="color:#808080;font-style:italic;"># Initialize original grid points
</span><span style="color:#808080;font-style:italic;"> </span>grid_dim = <span style="color:#0000ff;">10
</span><span style="color:#0000ff;"> </span>grid_spacing = <span style="color:#0000ff;">100
</span><span style="color:#0000ff;"> </span>num_points = grid_dim * grid_dim
points_original = np.zeros(<span style="color:#660099;">shape</span>=(num_points, <span style="color:#0000ff;">2</span>), <span style="color:#660099;">dtype</span>=np.float64)
<span style="color:#000080;font-weight:bold;">for </span>i <span style="color:#000080;font-weight:bold;">in </span>range(grid_dim):
<span style="color:#000080;font-weight:bold;">for </span>j <span style="color:#000080;font-weight:bold;">in </span>range(grid_dim):
points_original[i*grid_dim+j] = (i*grid_spacing, j*grid_spacing)
<span style="color:#808080;font-style:italic;"># Compute transformed grid points
</span><span style="color:#808080;font-style:italic;"> </span>affine_mat = np.array([[-<span style="color:#0000ff;">0.5</span>, -<span style="color:#0000ff;">0.86</span>, <span style="color:#0000ff;">100</span>], [<span style="color:#0000ff;">0.86</span>, -<span style="color:#0000ff;">0.5</span>, <span style="color:#0000ff;">100</span>]]) <span style="color:#808080;font-style:italic;"># createAffineMatrix(120, 1, 1, 100, 100)
</span><span style="color:#808080;font-style:italic;"> </span>points_transformed = np.matmul(affine_mat[:<span style="color:#0000ff;">2</span>,:<span style="color:#0000ff;">2</span>], points_original.T).T + affine_mat[:<span style="color:#0000ff;">2</span>,<span style="color:#0000ff;">2</span>]
<span style="color:#808080;font-style:italic;"># Initialize PETSc objects
</span><span style="color:#808080;font-style:italic;"> </span>num_unknown = points_transformed.size
mat_shape = (num_unknown, num_unknown)
A = PETSc.Mat()
A.createAIJ(<span style="color:#660099;">size</span>=mat_shape, <span style="color:#660099;">comm</span>=PETSc.COMM_WORLD)
A.setUp()
x, f = A.createVecs()
options = PETSc.Options()
options.setValue(<span style="color:#008000;font-weight:bold;">"-snes_qn_type"</span>, <span style="color:#008000;font-weight:bold;">"lbfgs"</span>) <span style="color:#808080;font-style:italic;"># broyden/lbfgs
</span><span style="color:#808080;font-style:italic;"> </span>options.setValue(<span style="color:#008000;font-weight:bold;">"-snes_qn_scale_type"</span>, <span style="color:#008000;font-weight:bold;">"none"</span>) <span style="color:#808080;font-style:italic;"># none, diagonal, scalar, jacobian,
</span><span style="color:#808080;font-style:italic;"> </span>options.setValue(<span style="color:#008000;font-weight:bold;">"-snes_monitor"</span>, <span style="color:#008000;font-weight:bold;">""</span>)
<span style="color:#808080;font-style:italic;"># options.setValue("-snes_view", "")
</span><span style="color:#808080;font-style:italic;"> </span>options.setValue(<span style="color:#008000;font-weight:bold;">"-snes_converged_reason"</span>, <span style="color:#008000;font-weight:bold;">""</span>)
options.setFromOptions()
snes = PETSc.SNES()
snes.create(PETSc.COMM_WORLD)
snes.setType(<span style="color:#008000;font-weight:bold;">"qn"</span>)<span style="color:#808080;font-style:italic;">
</span><span style="color:#808080;font-style:italic;"> </span>snes.setFunction(fill_function, f, <span style="color:#660099;">args</span>=(points_original,))
snes.setJacobian(fill_jacobian, A, <span style="color:#000080;">None</span>, <span style="color:#660099;">args</span>=(points_original,))
snes_pc = snes.getNPC() <span style="color:#808080;font-style:italic;"># Inner snes instance (newtonls by default!)
</span><span style="color:#808080;font-style:italic;"> # snes_pc.setType("ngmres")
</span><span style="color:#808080;font-style:italic;"> </span>snes.setFromOptions()
ksp = snes_pc.getKSP()
ksp.setType(<span style="color:#008000;font-weight:bold;">"cg"</span>)
ksp.setTolerances(<span style="color:#660099;">rtol</span>=<span style="color:#0000ff;">1e-10</span>, <span style="color:#660099;">max_it</span>=<span style="color:#0000ff;">40000</span>)
pc = ksp.getPC()
pc.setType(<span style="color:#008000;font-weight:bold;">"asm"</span>)
ksp.setFromOptions()
x.setArray(points_transformed.ravel())
snes.solve(<span style="color:#000080;">None</span>, x)</pre>
<br>
<font face="Arial" color="Black" size="1"><br>
------------------------------------------------------------------------------------------------<br>
------------------------------------------------------------------------------------------------<br>
Forschungszentrum Juelich GmbH<br>
52425 Juelich<br>
Sitz der Gesellschaft: Juelich<br>
Eingetragen im Handelsregister des Amtsgerichts Dueren Nr. HR B 3498<br>
Vorsitzender des Aufsichtsrats: MinDir Volker Rieke<br>
Geschaeftsfuehrung: Prof. Dr.-Ing. Wolfgang Marquardt (Vorsitzender),<br>
Karsten Beneke (stellv. Vorsitzender), Prof. Dr.-Ing. Harald Bolt<br>
------------------------------------------------------------------------------------------------<br>
------------------------------------------------------------------------------------------------<br>
<br>
</font>
</body>
</html>