[petsc-users] Singular jocabian using SNES
Marcel Huysegoms
m.huysegoms at fz-juelich.de
Wed Mar 24 05:19:13 CDT 2021
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.
However, I believe there is a *wrapper function missing in the petsc4py
library*.
When I perform (from the attached python script):
tao = PETSc.TAO()
tao.create(PETSc.COMM_WORLD)
tao.setType("brgn")
tao.setFromOptions()
tao.setResidual(fill_function, f,args=(points_original,))
tao.setJacobian(fill_jacobian, A,args=(points_original,))
tao.solve(x)
... I get the following error message
petsc4py.PETSc.Error: error code 58
[0] TaoSolve() line 215 in /opt/petsc/src/tao/interface/taosolver.c
[0] TaoSetUp() line 269 in /opt/petsc/src/tao/interface/taosolver.c
[0] TaoSetUp_BRGN() line 242 in
/opt/petsc/src/tao/leastsquares/impls/brgn/brgn.c
[0] Operation done in wrong order
[0] *TaoSetResidualJacobianRoutine() must be called before setup!*
The TAO documentation states, for solving LLSQ problems (i.e. using
BRGN) the two functions "setResidualRoutine()" and
"setJacobianResidualRoutine()" need to be called.
The first one is invoked by petsc4py.setResidual(). However there is no
petsc4py function that invokes the second routine.
Looking into the source code of petsc4py (TAO.pyx), there is only
setJacobian() which invokes ToaSetJacobianRoutine() but not
ToaSetJacobian*Residual*Routine()
https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Tao/TaoSetJacobianRoutine.html
https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Tao/TaoSetJacobianResidualRoutine.html
Am I missing something or is petsc4py actually lacking a wrapper
function like setJacobianResidual()?
Best regards,
Marcel
On 23.03.21 20:42, Matthew Knepley wrote:
> On Tue, Mar 23, 2021 at 12:39 PM Marcel Huysegoms
> <m.huysegoms at fz-juelich.de <mailto: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 sys
> import petsc4py
> import numpyas np
>
> petsc4py.init(sys.argv)
> from petsc4pyimport 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 iin 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 iin range(grid_dim):
> for jin 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/>
--
Marcel Huysegoms
Institut für Neurowissenschaften und Medizin (INM-1)
Forschungszentrum Jülich GmbH
52425 Jülich
Telefon: +49 2461 61 3678
Email: m.huysegoms at fz-juelich.de
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20210324/73622f4f/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: tao_test.py
Type: text/x-python
Size: 2096 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20210324/73622f4f/attachment.py>
More information about the petsc-users
mailing list