<div dir="ltr"><div dir="ltr">On Tue, Mar 23, 2021 at 12:39 PM Marcel Huysegoms <<a href="mailto:m.huysegoms@fz-juelich.de">m.huysegoms@fz-juelich.de</a>> wrote:<br></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<div>
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 href="https://scicomp.stackexchange.com/questions/21781/newtons" target="_blank">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></div></blockquote><div><br></div><div>SNES is only for solving systems of nonlinear equations. If you want optimization (least-square, etc.) then you want to formulate your</div><div>problem in the TAO interface. It has quasi-Newton methods for those problems, and other methods as well. That is where I would start.</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div>
Many thanks in advance,<br>
Marcel<br>
<br>
-----------------------------------------<br>
<br>
<pre style="background-color:rgb(241,241,234);color:rgb(0,0,0);font-family:"DejaVu Sans Mono""><span style="color:rgb(0,0,128);font-weight:bold">import </span>sys
<span style="color:rgb(0,0,128);font-weight:bold">import </span>petsc4py
<span style="color:rgb(0,0,128);font-weight:bold">import </span>numpy <span style="color:rgb(0,0,128);font-weight:bold">as </span>np
petsc4py.init(sys.argv)
<span style="color:rgb(0,0,128);font-weight:bold">from </span>petsc4py <span style="color:rgb(0,0,128);font-weight:bold">import </span>PETSc
<span style="color:rgb(0,0,128);font-weight:bold">def </span>fill_function(<span style="color:rgb(128,128,128)">snes</span>, x, f, points_original):
x_values = x.getArray(<span style="color:rgb(102,0,153)">readonly</span>=<span style="color:rgb(0,0,128)">True</span>)
diff_vectors = points_original.ravel() - x_values
f_values = np.square(diff_vectors)
<span style="color:rgb(128,128,128);font-style:italic"># f_values[10] = 0
</span><span style="color:rgb(128,128,128);font-style:italic"> </span>f.setValues(np.arange(f_values.size), f_values)
f.assemble()
<span style="color:rgb(0,0,128);font-weight:bold">def </span>fill_jacobian(<span style="color:rgb(128,128,128)">snes</span>, x, <span style="color:rgb(128,128,128)">J</span>, P, points_original):
x_values = x.getArray(<span style="color:rgb(102,0,153)">readonly</span>=<span style="color:rgb(0,0,128)">True</span>)
points_original_flat = points_original.ravel()
deriv_values = -<span style="color:rgb(0,0,255)">2</span>*(points_original_flat - x_values)
<span style="color:rgb(128,128,128);font-style:italic"># deriv_values[10] = 0
</span><span style="color:rgb(128,128,128);font-style:italic"> </span><span style="color:rgb(0,0,128);font-weight:bold">for </span>i <span style="color:rgb(0,0,128);font-weight:bold">in </span>range(x_values.size):
P.setValue(i, i, deriv_values[i])
<span style="color:rgb(128,128,128);font-style:italic"># print(deriv_values)
</span><span style="color:rgb(128,128,128);font-style:italic"> </span>P.assemble()
<span style="color:rgb(128,128,128);font-style:italic"># ---------------------------------------------------------------------------------------------
</span><span style="color:rgb(128,128,128);font-style:italic">
</span><span style="color:rgb(0,0,128);font-weight:bold">if </span>__name__ == <span style="color:rgb(0,128,0);font-weight:bold">'__main__'</span>:
<span style="color:rgb(128,128,128);font-style:italic"># Initialize original grid points
</span><span style="color:rgb(128,128,128);font-style:italic"> </span>grid_dim = <span style="color:rgb(0,0,255)">10
</span><span style="color:rgb(0,0,255)"> </span>grid_spacing = <span style="color:rgb(0,0,255)">100
</span><span style="color:rgb(0,0,255)"> </span>num_points = grid_dim * grid_dim
points_original = np.zeros(<span style="color:rgb(102,0,153)">shape</span>=(num_points, <span style="color:rgb(0,0,255)">2</span>), <span style="color:rgb(102,0,153)">dtype</span>=np.float64)
<span style="color:rgb(0,0,128);font-weight:bold">for </span>i <span style="color:rgb(0,0,128);font-weight:bold">in </span>range(grid_dim):
<span style="color:rgb(0,0,128);font-weight:bold">for </span>j <span style="color:rgb(0,0,128);font-weight:bold">in </span>range(grid_dim):
points_original[i*grid_dim+j] = (i*grid_spacing, j*grid_spacing)
<span style="color:rgb(128,128,128);font-style:italic"># Compute transformed grid points
</span><span style="color:rgb(128,128,128);font-style:italic"> </span>affine_mat = np.array([[-<span style="color:rgb(0,0,255)">0.5</span>, -<span style="color:rgb(0,0,255)">0.86</span>, <span style="color:rgb(0,0,255)">100</span>], [<span style="color:rgb(0,0,255)">0.86</span>, -<span style="color:rgb(0,0,255)">0.5</span>, <span style="color:rgb(0,0,255)">100</span>]]) <span style="color:rgb(128,128,128);font-style:italic"># createAffineMatrix(120, 1, 1, 100, 100)
</span><span style="color:rgb(128,128,128);font-style:italic"> </span>points_transformed = np.matmul(affine_mat[:<span style="color:rgb(0,0,255)">2</span>,:<span style="color:rgb(0,0,255)">2</span>], points_original.T).T + affine_mat[:<span style="color:rgb(0,0,255)">2</span>,<span style="color:rgb(0,0,255)">2</span>]
<span style="color:rgb(128,128,128);font-style:italic"># Initialize PETSc objects
</span><span style="color:rgb(128,128,128);font-style:italic"> </span>num_unknown = points_transformed.size
mat_shape = (num_unknown, num_unknown)
A = PETSc.Mat()
A.createAIJ(<span style="color:rgb(102,0,153)">size</span>=mat_shape, <span style="color:rgb(102,0,153)">comm</span>=PETSc.COMM_WORLD)
A.setUp()
x, f = A.createVecs()
options = PETSc.Options()
options.setValue(<span style="color:rgb(0,128,0);font-weight:bold">"-snes_qn_type"</span>, <span style="color:rgb(0,128,0);font-weight:bold">"lbfgs"</span>) <span style="color:rgb(128,128,128);font-style:italic"># broyden/lbfgs
</span><span style="color:rgb(128,128,128);font-style:italic"> </span>options.setValue(<span style="color:rgb(0,128,0);font-weight:bold">"-snes_qn_scale_type"</span>, <span style="color:rgb(0,128,0);font-weight:bold">"none"</span>) <span style="color:rgb(128,128,128);font-style:italic"># none, diagonal, scalar, jacobian,
</span><span style="color:rgb(128,128,128);font-style:italic"> </span>options.setValue(<span style="color:rgb(0,128,0);font-weight:bold">"-snes_monitor"</span>, <span style="color:rgb(0,128,0);font-weight:bold">""</span>)
<span style="color:rgb(128,128,128);font-style:italic"># options.setValue("-snes_view", "")
</span><span style="color:rgb(128,128,128);font-style:italic"> </span>options.setValue(<span style="color:rgb(0,128,0);font-weight:bold">"-snes_converged_reason"</span>, <span style="color:rgb(0,128,0);font-weight:bold">""</span>)
options.setFromOptions()
snes = PETSc.SNES()
snes.create(PETSc.COMM_WORLD)
snes.setType(<span style="color:rgb(0,128,0);font-weight:bold">"qn"</span>)<span style="color:rgb(128,128,128);font-style:italic">
</span><span style="color:rgb(128,128,128);font-style:italic"> </span>snes.setFunction(fill_function, f, <span style="color:rgb(102,0,153)">args</span>=(points_original,))
snes.setJacobian(fill_jacobian, A, <span style="color:rgb(0,0,128)">None</span>, <span style="color:rgb(102,0,153)">args</span>=(points_original,))
snes_pc = snes.getNPC() <span style="color:rgb(128,128,128);font-style:italic"># Inner snes instance (newtonls by default!)
</span><span style="color:rgb(128,128,128);font-style:italic"> # snes_pc.setType("ngmres")
</span><span style="color:rgb(128,128,128);font-style:italic"> </span>snes.setFromOptions()
ksp = snes_pc.getKSP()
ksp.setType(<span style="color:rgb(0,128,0);font-weight:bold">"cg"</span>)
ksp.setTolerances(<span style="color:rgb(102,0,153)">rtol</span>=<span style="color:rgb(0,0,255)">1e-10</span>, <span style="color:rgb(102,0,153)">max_it</span>=<span style="color:rgb(0,0,255)">40000</span>)
pc = ksp.getPC()
pc.setType(<span style="color:rgb(0,128,0);font-weight:bold">"asm"</span>)
ksp.setFromOptions()
x.setArray(points_transformed.ravel())
snes.solve(<span style="color:rgb(0,0,128)">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>
</div>
</blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div><div><br></div><div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>