<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
</head>
<body>
thanks a lot for your suggestion! The TAO interface seems to be
perfectly suited for the problem and I wasn't aware of its
existence.<br>
<br>
However, I believe there is a <b>wrapper function missing in the
petsc4py library</b>. <br>
When I perform (from the attached python script):<br>
<pre style="background-color:#f1f1ea;color:#000000;font-family:'DejaVu Sans Mono';font-size:10,5pt;">tao = PETSc.TAO()
tao.create(PETSc.COMM_WORLD)
tao.setType(<span style="color:#008000;font-weight:bold;">"brgn"</span>)
tao.setFromOptions()
tao.setResidual(fill_function, f, <span style="color:#660099;">args</span>=(points_original,))
tao.setJacobian(fill_jacobian, A, <span style="color:#660099;">args</span>=(points_original,))
tao.solve(x)
<span style="color:#808080;font-style:italic;"></span></pre>
<br>
... I get the following error message<br>
<br>
petsc4py.PETSc.Error: error code 58<br>
[0] TaoSolve() line 215 in /opt/petsc/src/tao/interface/taosolver.c<br>
[0] TaoSetUp() line 269 in /opt/petsc/src/tao/interface/taosolver.c<br>
[0] TaoSetUp_BRGN() line 242 in
/opt/petsc/src/tao/leastsquares/impls/brgn/brgn.c<br>
[0] Operation done in wrong order<br>
[0] <b>TaoSetResidualJacobianRoutine() must be called before setup!</b><br>
<br>
The TAO documentation states, for solving LLSQ problems (i.e. using
BRGN) the two functions "setResidualRoutine()" and
"setJacobianResidualRoutine()" need to be called.<br>
<br>
The first one is invoked by petsc4py.setResidual(). However there is
no petsc4py function that invokes the second routine.<br>
Looking into the source code of petsc4py (TAO.pyx), there is only
setJacobian() which invokes ToaSetJacobianRoutine() but not
ToaSetJacobian<b>Residual</b>Routine()<br>
<br>
<a class="moz-txt-link-freetext" href="https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Tao/TaoSetJacobianRoutine.html">https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Tao/TaoSetJacobianRoutine.html</a><br>
<a class="moz-txt-link-freetext" href="https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Tao/TaoSetJacobianResidualRoutine.html">https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Tao/TaoSetJacobianResidualRoutine.html</a><br>
<br>
Am I missing something or is petsc4py actually lacking a wrapper
function like setJacobianResidual()?<br>
<br>
Best regards,<br>
Marcel<br>
<br>
<span class="subject"></span><br>
<span class="subject"></span><br>
<br>
<br>
<div class="moz-cite-prefix">On 23.03.21 20:42, Matthew Knepley
wrote:<br>
</div>
<blockquote type="cite"
cite="mid:CAMYG4GnmkxX2A-HOn1fwigjgyTh919y+YH1N53-0aKJYumXR-Q@mail.gmail.com">
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<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"
moz-do-not-send="true">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" moz-do-not-send="true">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 size="1" face="Arial" color="Black"><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" moz-do-not-send="true">https://www.cse.buffalo.edu/~knepley/</a><br>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</blockquote>
<br>
<pre class="moz-signature" cols="72">--
Marcel Huysegoms
Institut für Neurowissenschaften und Medizin (INM-1)
Forschungszentrum Jülich GmbH
52425 Jülich
Telefon: +49 2461 61 3678
Email: <a class="moz-txt-link-abbreviated" href="mailto:m.huysegoms@fz-juelich.de">m.huysegoms@fz-juelich.de</a></pre>
</body>
</html>