I have an unconstrained minimization problem of many variables for which I know the gradient exactly. I turned to the conjugate gradient method contained in scipy.optimize.minimize (which uses the Polak-Ribiere algorithm), but it throws a LineSearchError when I try to converge the algorithm beyond the square root of machine precision. This seems to be a common occurence with a certain class of line-search algorithms.

The square root of machine precision is not enough for my purposes. Is there a robust algorithm available that uses an approximate line-search, or something similar, which does enable one to converge to machine precision?

  • 3
    $\begingroup$ What exactly do you mean with machine precision? Do you mean single/double precision? Why are you looking to converge to the square root of that? $\endgroup$ – MPIchael Sep 25 '19 at 14:26
  • 2
    $\begingroup$ How well conditioned is your problem? Do you know the condition number at the optimal solution? $\endgroup$ – Brian Borchers Sep 25 '19 at 16:50
  • $\begingroup$ I am looking to converge beyond sqrt(eps) (double precision) because I am after a scalar which contains two contractions with the variable vector at a local minimum, meaning the precision of that scalar will be again the square root of the precision of the variables. That turns out to be ~1e-4, and that's insufficient for my problem. @BrianBorchers I don't know the condition number, except that it should be fairly well behaved. $\endgroup$ – Kappie001 Sep 25 '19 at 17:09
  • 1
    $\begingroup$ If the condtion number is $10^{k}$, then the best you can hope for in double precision is $16-k$ digits of relative precision in the answer. $\endgroup$ – Brian Borchers Sep 25 '19 at 17:55
  • 2
    $\begingroup$ I have no reason to believe that there's anything wrong with the code that you're using. A good working hypothesis is that your problem has a condition number of around $10^8$ and that the scipy function is doing as well as any routine could using double precision. $\endgroup$ – Brian Borchers Sep 25 '19 at 18:57

I totally agree with the discussion in the comments: it is quite likely that your problem has a large enough condition number that you are seeing the problems converging beyond $\sqrt{\epsilon_\text{mach}}$.

Without knowing anything about the problem in particular, one computationally costly solution is to use higher (say, arbitrary) precision floating-point representation. While it might slow down the performance significantly, it should be relatively easy to implement.

Since you are working with Python, I suggest looking into mpmath. This is a Python library offering arbitrary precision. In particular, you should look into their implementation of root finding and optimization:

class mpmath.calculus.optimization.MDNewton(ctx, f, x0, **kwargs`)

Note, they still recommend obtaining the initial guess from SciPy (to improve performance), and in your case, you might use mpmath simply to refine the solution obtained by SciPy. Notice, that you probably will have to store most of the things in the arbitrary precision.

| cite | improve this answer | |
  • $\begingroup$ No, the matrix definitely does not have a condition number of 10^{8}. I ended up writing this code to Python, and it converges to machine precision. It uses approximate line search. $\endgroup$ – Kappie001 Sep 26 '19 at 15:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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