[petsc-users] Fwd: Implementing a homotopy solver
Matthew Knepley
knepley at gmail.com
Fri Oct 5 15:42:00 CDT 2018
On Fri, Oct 5, 2018 at 4:15 PM zakaryah <zakaryah at gmail.com> wrote:
> 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?
>
I do not claim to fully understand the above (it would be easier for me to
write the equations explicitly on their own lines). However, the issue with
nullspaces is soluble I think. So you say that
J \in R^{m \cross (m+1)} is full rank
so that its nullspace has dimension 1 and is spanned by u,
J u = 0.
Now you augment J to get a square matrix A by adding a row,
/ J \ u = / J u \ = / 0 \
\ v^T / \ v^T u / \ 0 /
if v is orthogonal to u. Thus, just choose v from the m-1 dimensional space
orthogonal to u.
Thanks,
Matt
> 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/>
>
>
>
--
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/aa0807d4/attachment-0001.html>
More information about the petsc-users
mailing list