5
$\begingroup$

I am using the Newton's method to solve $3\times3$ systems.

For some particular cases, it turns out that at a given iteration, the Jacobian matrix cannot be inverted and that its determinant is very close to zero (looking at the matrix, there are terms that are around 1e+0 and others that are 1e-15).

After investigations, it is clear that one variable has no influence on the system when it is close to the solution.

What is the most clever thing to do with such an issue ?

I would like to have an algorithm that can adapt itself to such situations when they happen.

EDIT

It is about an optical optimisation problem. The point is to add a surface to an optical system so that it fits some optical properties. The Newton's method finds the roots of a function that takes as input the parameters of the surface and outputs the differences between the optical properties computed and the targeted properties.

I noticed that if the system is complex enough, then we have convergence. But if the system is too simple, the surface is more spherical and the Jacobian goes to very small values, because some parameters like astigmatism axis become influence-less.

Improve this question
$\endgroup$
7
  • 4
    $\begingroup$ If these systems arise from an optimization problem, you could look at standard globalization approaches such as line searches or trust region methods. However, Newton's method just doesn't work if the derivative vanishes in a neighborhood of the solution. You could then try a quasi-Newton method such as Broyden's (or BFGS). $\endgroup$
    – Christian Clason
    Jan 11, 2016 at 11:20
  • $\begingroup$ Related question (possibly duplicate?): scicomp.stackexchange.com/questions/4781/… $\endgroup$
    – Christian Clason
    Jan 11, 2016 at 11:41
  • 1
    $\begingroup$ It seams that I have found a simple trick : when the Jacobian determinant is too small, I arbitrary replace the column that causes the problem by a $(0,...,0,1)$ column. So it gives a fake "action on the system" to the variable. Because the variable is influence-less on the system, whatever the solution obtained, the result fits the equations.... $\endgroup$
    – Sylvain B.
    Jan 11, 2016 at 14:12
  • $\begingroup$ Thanks Christian for showing interests in my question. I do not understand clearly your last comment. Could you reformulate please ? Could you help me to build the proof? $\endgroup$
    – Sylvain B.
    Jan 18, 2016 at 12:25
  • $\begingroup$ I guess it is more a root finding than a minimization, because the solution may always exist. The Newton's method works with an algorithm that computes the optical properties. The Jacobian is evaluated at each iteration by finite differences using the algorithm. The three parameters in the loop describes the surface curvature properties (principal axes and values). And the algorithm crashes when the surface is too close to a sphere. Then, the influence less variable is the axis (projected on a plane so it is a angle). It is the column associated to this variable that I replace arbitrary. $\endgroup$
    – Sylvain B.
    Jan 18, 2016 at 13:17

1 Answer 1

Reset to default
3
$\begingroup$

From what you describe, you have an ill-posed problem: The solution is not unique (not even locally). A standard way of dealing with this is the following: Instead of trying to solve $F(x)=y$, where $x$ is your vector of geometrical parameters, $y$ the vector of optical parameters, and $F$ the mapping that computes the latter via the former, you minimize $$ J(x) = \frac12\|F(x)-y\|^2 + \frac\alpha2\|x\|^2$$ for some (small) $\alpha>0$. The last term ensures (local) uniqueness: Among all solutions of $F(x)=y$, it will pick the one with minimal norm (in your case, if any angle will give the same surface, the minimizer will have zero angle).

You can then compute a minimizer using any optimization method (e.g., BFGS with line search globalization), as described in the book by Nocedal and Wright Numerical Optimization.

For most methods, you need the gradient of $J(x)$, which is given by $$\nabla J(x) = \nabla F(x)^T (F(x)-y) + \alpha x,$$ where $\nabla F(x)^T$ is the transpose of the Jacobian. If it is at all possible, you should try to compute the Jacobian analytically and (approximately) evaluate that. (If you add the mathematical description of your $F$, we might be able to help with that.)

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thank you for this answer. I would even more use $\alpha = 0$ because i don't need a "unique" solution. I could check that the algorithm is robust on many test cases. I never mentioned performance issues in my question, but the $F$ function is a time-consuming algorithm that makes this optimization method $10\times$slower than a simple root finding algorithm using Broyden's method, but I will remember this variant trick. PS : sorry for answering so late, I was quite busy and I couldn't find time to workout your suggestion. $\endgroup$
    – Sylvain B.
    Mar 11, 2016 at 16:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.