[petsc-users] Fwd: Implementing a homotopy solver

Fri Oct 5 15:14:35 CDT 2018

Thanks Matt - I was trying to write it in linear algebra notation, but the
nature of the problem might produce some confusion.  I was also originally
a bit confused by a fairly major typo in the paper I was following, but I'm
almost sure I have it clear now.  I'll try to be concise - the meat is in
the last two paragraphs.

My system of nonlinear equations, F(x)=0, has dimension m.  Traditional
methods for solving, like Newton's method, tend to reach points at which
the Jacobian becomes singular.  I'm trying the homotopy trick, which
involves introducing an additional scalar variable, s, and solving the
system H(x,s) = sF(x) + (1-s)(x-x_0), starting at s=0 and following the
zero curve to, hopefully, a solution at s=1.

The homotopy map H only has m components, so its Jacobian J is m x (m+1),
and has rank m.  I need to find a vector u in its null space: J u = 0.
Note that u is size m+1.  Rather than solve J u = 0, I make an augmented
matrix A, which is (m+1)x(m+1), and consists of J in its first m rows, then
a somewhat arbitrary vector in the last row.  Then the vector u is the
solution to A u = b, where b is zero everywhere except its last element.
This solve seems to be working fine.

The SNES should then solve for a particular solution y of the Newton
update, J y = - H(y).  A (unique) particular solution can also be found by
solving A y = [-H(y) 0]^T.  Finally - and this is the step I can't figure
out how to implement properly, I want to find the norm-minimizing solution
z for the Newton update:
z = y - u(y^T u)/(u^T u).

I doubt there is any way to do this elegantly within the existing routines,
because u is not in the nullspace of A.  Rather, it's in the nullspace of
J, which is the submatrix I care about.  I think what I need is a way to
just apply the orthogonal projection in the preceding paragraph, after the
SNES solves the linear system but before the Newton update is applied.  I
don't know the SNES internals well at all.  Is it possible to hook this in
somehow, or put it in an existing function?

Thanks again for all your help.

On Oct 1, 2018 6:49 AM, "Matthew Knepley" <knepley at gmail.com> wrote:

On Sun, Sep 30, 2018 at 6:06 PM zakaryah <zakaryah at gmail.com> wrote:

> OK, thanks.
>
> I'm using a composite DM, DM_packer.  To make a separate KSP, I did the
> following in the FormJacobian() routine, after assembling the Jacobian
> matrix A and the RHS for the tangent vector, b:
>
>    - KSPCreate(PETSC_COMM_WORLD,&ksp)
>    - KSPSetDM(ksp,DM_packer)
>    - KSPSetDMActive(ksp,PETSC_FALSE) - because I want to set the operators
>    - KSPSetOperators(ksp,A,P)
>    - VecSet(n,0) - set initial guess to zero
>    - KSPSolve(ksp,b,n) - this solve works correctly
>    - VecNormalize(n,NULL)
>    -
>    MatNullSpaceCreate(PetscObjectComm((PetscObject)A),PETSC_FALSE,1,&n,&nullsp)
>    - MatSetNullSpace(A,nullsp)
>
> Then, with -snes_type newtonls -pc_type none -ksp_monitor
> -ksp_monitor_true_residual -ksp_view, the output for the KSPSolve described
> above, i.e. for the tangent vector, looks correct:
>
>   0 KSP preconditioned resid norm 1.000000000000e+03 true resid norm
> 1.000000000000e+03 ||r(i)||/||b|| 1.000000000000e+00
> ...
>
> 185 KSP preconditioned resid norm 9.900713131874e-03 true resid norm
> 9.900713131904e-03 ||r(i)||/||b|| 9.900713131904e-06
>
> Linear solve converged due to CONVERGED_RTOL iterations 185
>
> KSP Object: 1 MPI processes
>
>   type: gmres
>
>     restart=30, using Classical (unmodified) Gram-Schmidt
> Orthogonalization with no iterative refinement
>
>     happy breakdown tolerance 1e-30
>
>   maximum iterations=10000, initial guess is zero
>
>   tolerances:  relative=1e-05, absolute=1e-50, divergence=10000.
>
>   left preconditioning
>
>   using PRECONDITIONED norm type for convergence test
>
> PC Object: 1 MPI processes
>
>   type: none
>
>   linear system matrix = precond matrix:
>
>   Mat Object: 1 MPI processes
>
>     type: seqaij
>
>     rows=78247, cols=78247
>
>     total: nonzeros=6063481, allocated nonzeros=6063481
>
>     total number of mallocs used during MatSetValues calls =0
>
>       using I-node routines: found 26083 nodes, limit used is 5
>
>
> But the next KSP, i.e. the one in the SNES, doesn't converge in the true
> residual:
>
>
>   0 KSP preconditioned resid norm 2.045599092896e-04 true resid norm
> 2.803828296212e-04 ||r(i)||/||b|| 1.000000000000e+00
>
> 191 KSP preconditioned resid norm 2.009941278534e-09 true resid norm
> 1.636010142734e-04 ||r(i)||/||b|| 5.834915586465e-01
>
>
> KSP Object: 1 MPI processes
>
>   type: gmres
>
>     restart=30, using Classical (unmodified) Gram-Schmidt
> Orthogonalization with no iterative refinement
>
>     happy breakdown tolerance 1e-30
>
>   maximum iterations=10000, initial guess is zero
>
>   tolerances:  relative=1e-05, absolute=1e-50, divergence=10000.
>
>   left preconditioning
>
>   using PRECONDITIONED norm type for convergence test
>
> PC Object: 1 MPI processes
>
>   type: none
>
>   linear system matrix = precond matrix:
>
>   Mat Object: 1 MPI processes
>
>     type: seqaij
>
>     rows=78247, cols=78247
>
>     total: nonzeros=6063481, allocated nonzeros=6063481
>
>     total number of mallocs used during MatSetValues calls =0
>
>       has attached null space
>
>       using I-node routines: found 26083 nodes, limit used is 5
>
>
> My guess about what's going on is that the tangent vector n isn't really
> in the nullspace of A.
>

I do not understand the rest, but this is the problem. Maybe it would make
more sense if you wrote things
in linear algebraic notation.

  Thanks,

     Matt

> Rather, it's in the nullspace of the m x (m+1) submatrix of A.  So, An=c
> e_{m+1}, where c is an arbitrary constant and e_{m+1} is the m+1 th basis
> vector.  The nonlinear function also has an added row, F_{m+1}, which is
> set to zero in the FormFunction() routine.  I don't care about the value of
> F_{m+1}, but I suppose that if it's included in the true residual, and NOT
> in the "preconditioned" residual, even with pc_type none, then it will be
> tricky to diagnose the performance.  Should I be using my own routine to
> evaluate the residual, so that F_{m+1} is not included?
>

-- 
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/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20181005/59e2049d/attachment-0001.html>