Issue 1:
http://docs.scipy.org/doc/scipy/reference/generated/scipy.cluster.hierarchy.linkage.html

The first couple of paragraphs say that y must be a condensed distance matrix, but the Parameters section says that y is optionally a data array. In practice the function does different things depending on the shape of y: if it appears 1d then it is treated as a condensed distance matrix, otherwise if it is 2d then it is treated as a data matrix (not as a compressed or uncompressed distance matrix). If it is not 1d or 2d then it falls off the if/elif and tries to return something that doesn't exist.

Here's the MATLAB function that hierarchy.linkage appears to try to clone:

http://www.mathworks.com/help/stats/linkage.html

Issue 2:
This submodule has broken my trust.

I am genuinely in awe of how it was included in scipy 1.0. I am tagging two heavy hitters to attempt to bring discussion to this issue since, in general, scipy.optimize issues seem to go without discussion or fixes. @pv @teoliphant

As I understand, an effort was made in the last few version of scipy leading up to and including 1.0 to unify the interface of the different optimization routines and eventually deprecate fmin_xxxx.

A fairly large issue with this is that many of the inputs and outputs in the scipy documentation for these functions are simply wrong.

Take, for example the BFGS method.

the documentation for disp is incorrect; convergence messages are not printed, only a blurb at the end. return_all appears in the call signature at the top and is undocumented. I looked to fmin_bfgs hoping its more specific documentation would be superior. The documentation for that is even more misleading. I hoped that the specified type for the constituents of allvecs being OptimizationResult would mean they hold cost function values and input vector values. They are simply input vector values. The documentation here is also misleading.

I also understand that many of the algorithms are wrapped code in another language, e.g. L-BFGS-B is written in FORTRAN77, and this may make their modification difficult. Earlier today I opened an issue that L-BFGS-B does not respect the maxiter input, #8372. There are other issues, e.g. lack of respect for user's input dtype because of hard coded doubles in the fortran code, for example.

However, the docs being misleading is a great travesty for the users of this module, which I believe is among the more important/widely used of scipy. I would work to improve them myself, but I do not feel that my knowledge of the inner workings of the module is great enough to do so without myself making mistakes.

Returning to L-BFGS-B, the wrapped fortran code appears to contain a variety of errors or behaviors that have a varied level of harmfulness. For example, this is an example output:

RUNNING THE L-BFGS-B CODE

* * *

Machine precision = 2.220D-16
N = 4 M = 10

At X0 0 variables are exactly at the bounds

At iterate 0 f= 2.75370D+01 |proj g|= 1.06893D+02

At iterate 1 f= 2.27149D+01 |proj g|= 1.76623D+02

At iterate 2 f= 1.20196D+01 |proj g|= 1.20733D+02

At iterate 3 f= 8.05081D+00 |proj g|= 8.74848D+01

At iterate 4 f= 5.00344D+00 |proj g|= 8.42253D+01

At iterate 5 f= 3.70940D+00 |proj g|= 1.01991D+02

At iterate 6 f= 1.57818D+00 |proj g|= 7.10948D+01

At iterate 7 f= 4.81956D-01 |proj g|= 3.76945D+01

At iterate 8 f= 6.81301D-02 |proj g|= 2.20645D+01

At iterate 9 f= 1.34639D-02 |proj g|= 8.79078D+00

* * *

Tit = total number of iterations
Tnf = total number of function evaluations
Tnint = total number of segments explored during Cauchy searches
Skip = number of BFGS updates skipped
Nact = number of active bounds at final generalized Cauchy point
Projg = norm of the final projected gradient
F = final function value

* * *

N Tit Tnf Tnint Skip Nact Projg F
4 9 12 1 0 0 8.791D+00 1.346D-02
F = 1.346394518877836E-002

CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH

Cauchy time 0.000E+00 seconds.
Subspace minimization time 0.000E+00 seconds.
Line search time 0.000E+00 seconds.

Total User time 0.000E+00 seconds.

the code writes "iterate" instead of "iteration" (a typo). This is an annoyance as it is not very professional, but it does not break anything. The separate timings for Cauchy, subspace minimization, and line search are also always zero. Same for user time. It would be nice for these to not be a filler values since this is misleading, but same thing - this is creature comfort problem, not an accuracy problem.

My application is in a topic known as wavefront sensing, or phase retrieval. We put in some number of parameters, and the landscape of the cost function is highly evolutionary as the optimizer traverses the space. As a result, when designing a cost function or debugging code it is desirable to know the state of the optimizer, even without knowing its internal mechanisms, as it walks the space. To do this I do the following:

Save this output text from L-BFGS-B, and scrape the f= values from it.

Use a callback that stores parameter vectors.

I expect these two things to have the same length always, and depend on cost[n] being associated with all_x[n] without registration errors between these two blocks of data.

In my particular application (I can share code on request, but 1000-2000 lines are called to implement my routine split across a few modules; there is no easy "def foo()" analog that I have found to reproduce this) I am having trouble understanding why my cost function appears to not be correlated, based on the output of l-bfgs-b, to the closeness of the current iteration input to truth:

This has cost me dozens of hours of debugging time, only to learn that something is wrong with L-BFGS-B that causes the output of options={'disp': True} to not be the the cost function at a given parameter vector.

Because I capture stdout, I am not able to use print. I have embedded a line in my cost function that appends the current parameter vector and cost to a .txt file and resorted to matching values in that to L-BFGS-B print output to see what is wrong. The text from that is below, with numbered lines

00. params = [ 0. 0. 0. 0.]	f(x) = 27.5369873872

01. params = [ 1.00000000e-08 0.00000000e+00 0.00000000e+00 0.00000000e+00]	f(x) = 27.5369884562

02. params = [ 0.00000000e+00 1.00000000e-08 0.00000000e+00 0.00000000e+00]	f(x) = 27.5369867324

03. params = [ 0.00000000e+00 0.00000000e+00 1.00000000e-08 0.00000000e+00]	f(x) = 27.5369870096

04. params = [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e-08]	f(x) = 27.5369875808

05. params = [-0.80769268 0.49476373 0.28537275 -0.14630029]	f(x) = 86.4418490547

06. params = [-0.80769267 0.49476373 0.28537275 -0.14630029]	f(x) = 86.4418490135

07. params = [-0.80769268 0.49476374 0.28537275 -0.14630029]	f(x) = 86.441849366

08. params = [-0.80769268 0.49476373 0.28537276 -0.14630029]	f(x) = 86.4418499781

09. params = [-0.80769268 0.49476373 0.28537275 -0.14630028]	f(x) = 86.4418482602

10. params = [-0.16406864 0.1005026 0.05796848 -0.02971834]	f(x) = 22.714940612

11. params = [-0.16406863 0.1005026 0.05796848 -0.02971834]	f(x) = 22.7149388457

12. params = [-0.16406864 0.10050261 0.05796848 -0.02971834]	f(x) = 22.7149405671

13. params = [-0.16406864 0.1005026 0.05796849 -0.02971834]	f(x) = 22.7149404753

14. params = [-0.16406864 0.1005026 0.05796848 -0.02971833]	f(x) = 22.7149421094

15. params = [-0.10087717 0.09158391 0.05575493 -0.06632038]	f(x) = 12.0195875306

16. params = [-0.10087716 0.09158391 0.05575493 -0.06632038]	f(x) = 12.0195863232

17. params = [-0.10087717 0.09158392 0.05575493 -0.06632038]	f(x) = 12.0195878653

18. params = [-0.10087717 0.09158391 0.05575494 -0.06632038]	f(x) = 12.0195883001

19. params = [-0.10087717 0.09158391 0.05575493 -0.06632037]	f(x) = 12.0195867159

20. params = [ 0.15842029 -0.01548538 -0.01632799 -0.07159999]	f(x) = 48.7174210299

21. params = [ 0.1584203 -0.01548538 -0.01632799 -0.07159999]	f(x) = 48.7174233933

22. params = [ 0.15842029 -0.01548537 -0.01632799 -0.07159999]	f(x) = 48.717421743

23. params = [ 0.15842029 -0.01548538 -0.01632798 -0.07159999]	f(x) = 48.7174193993

24. params = [ 0.15842029 -0.01548538 -0.01632799 -0.07159998]	f(x) = 48.7174245544

25. params = [-0.05512618 0.07269238 0.04303647 -0.06725192]	f(x) = 8.05081234605

26. params = [-0.05512617 0.07269238 0.04303647 -0.06725192]	f(x) = 8.05081192026

27. params = [-0.05512618 0.07269239 0.04303647 -0.06725192]	f(x) = 8.0508114712

28. params = [-0.05512618 0.07269238 0.04303648 -0.06725192]	f(x) = 8.05081281076

29. params = [-0.05512618 0.07269238 0.04303647 -0.06725191]	f(x) = 8.05081206742

30. params = [-0.02149653 0.09172596 0.0267816 -0.06652566]	f(x) = 5.0034401551

31. params = [-0.02149652 0.09172596 0.0267816 -0.06652566]	f(x) = 5.00344002938

32. params = [-0.02149653 0.09172597 0.0267816 -0.06652566]	f(x) = 5.00343999403

33. params = [-0.02149653 0.09172596 0.02678161 -0.06652566]	f(x) = 5.00344079846

34. params = [-0.02149653 0.09172596 0.0267816 -0.06652565]	f(x) = 5.00343931285

35. params = [ 0.03012534 0.12108955 -0.01347029 -0.04719621]	f(x) = 3.70939759272

36. params = [ 0.03012535 0.12108955 -0.01347029 -0.04719621]	f(x) = 3.70939819132

37. params = [ 0.03012534 0.12108956 -0.01347029 -0.04719621]	f(x) = 3.70939842652

38. params = [ 0.03012534 0.12108955 -0.01347028 -0.04719621]	f(x) = 3.70939740919

39. params = [ 0.03012534 0.12108955 -0.01347029 -0.0471962 ]	f(x) = 3.70939657281

40. params = [ 0.00300894 0.11308784 0.00337076 -0.03858502]	f(x) = 1.57818376805

41. params = [ 0.00300895 0.11308784 0.00337076 -0.03858502]	f(x) = 1.578184007

42. params = [ 0.00300894 0.11308785 0.00337076 -0.03858502]	f(x) = 1.57818379941

43. params = [ 0.00300894 0.11308784 0.00337077 -0.03858502]	f(x) = 1.57818414645

44. params = [ 0.00300894 0.11308784 0.00337076 -0.03858501]	f(x) = 1.5781830571

45. params = [-0.01561316 0.12143126 0.00100096 -0.01330919]	f(x) = 0.481955650202

46. params = [-0.01561315 0.12143126 0.00100096 -0.01330919]	f(x) = 0.481955273258

47. params = [-0.01561316 0.12143127 0.00100096 -0.01330919]	f(x) = 0.481955840221

48. params = [-0.01561316 0.12143126 0.00100097 -0.01330919]	f(x) = 0.481955839908

49. params = [-0.01561316 0.12143126 0.00100096 -0.01330918]	f(x) = 0.481955373719

50. params = [-0.00041042 0.12401147 -0.00485423 0.00207093]	f(x) = 0.0681301442419

51. params = [-0.00041041 0.12401147 -0.00485423 0.00207093]	f(x) = 0.0681301690496

52. params = [-0.00041042 0.12401148 -0.00485423 0.00207093]	f(x) = 0.0681301089522

53. params = [-0.00041042 0.12401147 -0.00485422 0.00207093]	f(x) = 0.0681299235971

54. params = [-0.00041042 0.12401147 -0.00485423 0.00207094]	f(x) = 0.0681302701249

55. params = [-0.00144057 0.12278377 -0.00124594 -0.00264166]	f(x) = 0.0134639451888

56. params = [-0.00144056 0.12278377 -0.00124594 -0.00264166]	f(x) = 0.0134639411113

57. params = [-0.00144057 0.12278378 -0.00124594 -0.00264166]	f(x) = 0.0134638572809

58. params = [-0.00144057 0.12278377 -0.00124593 -0.00264166]	f(x) = 0.0134639506651

59. params = [-0.00144057 0.12278377 -0.00124594 -0.00264165]	f(x) = 0.0134639172132

The callback creates the following list:

[array([ 0., 0., 0., 0.]),
array([-0.00144057, 0.12278377, -0.00124594, -0.00264166]),
array([-0.00144057, 0.12278377, -0.00124594, -0.00264166]),
array([-0.00144057, 0.12278377, -0.00124594, -0.00264166]),
array([-0.00144057, 0.12278377, -0.00124594, -0.00264166]),
array([-0.00144057, 0.12278377, -0.00124594, -0.00264166]),
array([-0.00144057, 0.12278377, -0.00124594, -0.00264166]),
array([-0.00144057, 0.12278377, -0.00124594, -0.00264166]),
array([-0.00144057, 0.12278377, -0.00124594, -0.00264166]),
array([-0.00144057, 0.12278377, -0.00124594, -0.00264166])]

I have reproduced this for completeness, but the important lines are 00, the initial guess, which matches iterate 0 in the output text from the optimizer. Iteration 1 matches any of line 12-14. The parameter vectors are similar to:

[-0.16406864 0.10050261 0.05796848 -0.02971834]

The callback contains this input array at that iteration:

[-0.00144057, 0.12278377, -0.00124594, -0.00264166]

Which would be associated with a cost function 0.013, not 22 and change. As a result, I am unable to trust the tools I have to diagnose a problem with them. The output of the L-BFGS-B print simply is not the value of the cost function for the parameter vector at that iteration. I believe the problem is something like as follows:

The callback is called at the end of a line search before beginning a new gradient hunt

The f= is written at the start of a line search

But I cannot trust the documentation of scipy to be correct that the callback is called at the end of an iteration (~= end of line search in BFGS), and I cannot trust lbfgsb.f to be error free in this way.

Are there plans to improve the optimization module in a "months not years" timespan? Is anyone willing to rewrite lbfgsb.f or even port it to C/C++ or pure python to solve these problems? In commercialized applications, it is often desirable for a user to keep their finger on the pulse by monitoring convergence and potentially stopping optimization, hand tuning the inputs, and then continuing the work. The architecture of optimize.minimize is prohibitive to this. It would be nice to see a context manager based implementation that allows interruption of the optimizer without losing the current state, as this sort of behavior is only possible at the moment through use of a complicated callback function.