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
<[email protected] <mailto:[email protected]>> 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: [email protected]
import sys
import petsc4py
import numpy as np
petsc4py.init(sys.argv)
from petsc4py import PETSc
def fill_function(tao, 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(tao, x, J, 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):
J.setValue(i, i, deriv_values[i])
# print(deriv_values)
J.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("-tao_monitor", "")
options.setValue("-tao_view", "")
options.setFromOptions()
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,))
x.setArray(points_transformed.ravel())
tao.setInitial(x)
tao.solve(x)
