[GSoC] ENH: New least-squares algorithms #5044

nmayorov · 2015-07-13T21:26:37Z

Hi! This is an extension of my previous PR #5019.

I moved some commits to the previous PR, now the first new commit is 2710986

It adds special adjustments to algorithms which allow to solve large problems, when Jacobian matrix is sparse (only matrix-vector products are computed). Also sparse Jacobians can be estimated fast by finite differencing, when each row contains only few non-zero elements.

UPDATE:

I think that's about it what I wanted to implement in least_squares function. I would appreciate more active feedback and code review.

Here are two examples of usage:

http://nbviewer.ipython.org/gist/nmayorov/6098a514cc3277d72dd7 - sparse features.
http://nbviewer.ipython.org/gist/nmayorov/dac97f3ed9d638043191 - robust loss functions.

I see that conceptually an example with interesting bounds usage is missing. Your suggestions in this directions are most welcome.

ev-br · 2015-07-22T18:01:10Z

scipy/optimize/least_squares.py

+                J = J.toarray()
+            elif not issparse(J):
+                J = np.atleast_2d(J)
+


So, this checks sparsity and warns on each iteration. Better move the checks, warnings and construction of an appropriate Jacobian out of jac_wrapped. Same for the else branch below.

The warning will appear only once. I thought about it and this approach seemed as the most clear for me. That's basically the reason to introduce jac_wrapped to handle all trivia in it.

If the Jacobian is know to be dense, why checking if issparse(J) repeatedly.

Because it does no harm and takes the least amount of code (most readable) :). But, as I suggested above, let's just remove these conversions and checks all together (if a user wants to convert to dense - it's trivial to do so in his code).

ev-br · 2015-07-22T18:05:28Z

scipy/optimize/least_squares.py

+    n = x0.size
+
+    if diff_step is None:
+        epsfcn = np.finfo(float).eps


Nitpick: The call to finfo(float) is repeated at least four times, twice at module level. Better do it once, store the result, and import.

Sure, better import it from ._lsq_common.

ev-br · 2015-07-22T18:14:45Z

scipy/optimize/least_squares.py

+    if jac in ['2-point', '3-point']:
+        if jac_sparsity is not None:
+            if method == 'lm':
+                warn("`jac_sparsity` is ignored for method='lm', dense "


I'd make this an error.

Motivation? To avoid heavy computations? I think you might be right, better to stop here.

On the other hand there are a lot of similar places: converting sparse Jacobian to dense is also looking for troubles. Maybe leave it to a user to choose correct options and only warn if something happening.

LM and sparsity do not play together, this can never be a useful thing.
Generally, with this many mysterious switches I think it makes sense to weed out obviously wrong things.

Yes, I think you are right! Let's also disallow sparse->dense conversion, i. e. be more strict and explicit. Agree?

I'm interested to know your opinion.

ev-br · 2015-07-22T20:16:19Z

scipy/optimize/_group_columns.pyx

+    return groups
+
+
+def group_sparse(int m, int n, int [:] indices, int [:] indptr):


Is using int for indices overflow-safe here? I'd unlikely someone is using this large matrices, but still.

Potentially might be not, but np.int32 is the type for sparse matrices indices and indptr, so it can be changed only on much more fundamental level.

I think this link is again relevant https://github.com/scikit-learn/scikit-learn/wiki/C-integer-types:-the-missing-manual

You should be careful about this. As of Scipy 0.14, sparse matrices use int32 unless the number of nonzero elements exceeds 2^31, at which point int64 is used. So this is probably fine given the application here, but it is possible to get 64 bits ints as indices from a sparse matrix.

Both indices and indptr are affected by that? Should I then cast the arrays to int64?

On the other hand, arrays with more than 2^31 elements are already too big for the problem.

I think the clean solution would be to use Cython's fused types here, but it looks like slicing a huge matrix to obtain a small one brings indices.dtype back to int32, so small matrices can be expected not to have int64 indices.

ev-br · 2015-07-22T20:53:22Z

Two random comments:

a link to https://nmayorov.wordpress.com could be added to some comments or a module docstring :-)
least_squares.py should be renamed to _least_squares.py

nmayorov · 2015-07-23T01:31:58Z

I've done some refactoring of Jacobian validation code. I think it's better now.

ev-br · 2015-07-23T15:02:45Z

Running python runtests -s optimize --coverage shows several untested branches: the quartic equation in solve_trust_region_2d, constrained Cauchy step in dogbox, and a couple of corner cases in trf's handling of the reflected step. While 100% coverage is not a goal, it's better to avoid leaving large holes in at least smoke testing.

ev-br · 2015-07-23T15:13:39Z

Re validation/wrapping of Jacobian: I still prefer to only do this validation once: nmayorov#1

charris · 2015-09-09T16:21:34Z

scipy/optimize/_lsq/common.py

+
+def solve_lsq_trust_region(n, m, uf, s, V, Delta, initial_alpha=None,
+                           rtol=0.01, max_iter=10):
+    """Solve a trust-region problem arising in least-squares minimization by


The summary lines should be kept < 80 characters total. So for this, maybe something like

"""Use MINPACK method to solve a trust-region least squares problem blah, blah """ Or some such, depending on what the function actually does, which is not clear; "problem" could refer to almost anything.

Where does "< 80 characters total" come from? I can find several methods in scipy which have multiline "Short summary" and it looks really well on the site. And this functions are not public actually.

It is said "trust-region problem", which is quite well known and specific in optimization.

it is inherited from the Python docstring convention https://www.python.org/dev/peps/pep-0257/, in combination with the pep 8 line length limitation.

It is hard to write short summaries, but the multiline ones are definitely a style violation. Travis doesn't check that though.

Agree that we should respect PEP Style recommendations (maybe not always completely). I will try to figure out one-liners.

nmayorov · 2015-09-13T18:52:04Z

I want to rename scaling='jac' to scaling='auto'. @ev-br like/dislike?

About the last example: maybe keep it? Yes doctring is long, but you don't need to read it all (if you don't want to). On the other hand, we have now examples covering all features. Matlab docpages are very long too.

My general plan:

Merge this PR.
Rebase/finish linear LSQ PR.
Work on tutorial / examples.
Finish curve_fit RP.

ev-br · 2015-09-14T21:27:56Z

Your plan sounds good :-).

On scaling='jac' or 'auto' I am ambivalent. 'jac' sounds a bit more specific, but that's not something I'm going to have a strong opinion about.

Re last example: it does feel a little out of place. It reads a little more narrative --- which is good!, but more appropriate for a tutorial. So maybe move it over there or #5233. (we need to advertise tutorials more!)

One little knot to tie is #5044 (comment). It's likely not a serious issue, but it'd be good to not leave holes where reasonable. I've to admit I've no idea. Maybe @ewmoore or @larsmans can weigh in on this?

nmayorov · 2015-09-15T20:59:29Z

About huge sparse matrices: I just couldn't test in on my machine, as I said, it is unlikely to create troubles (time / our of ram issues will come first).

Updated: I meant I couldn't test in on my machine. Sorry about that.

larsmans · 2015-09-15T22:03:44Z

scipy/optimize/_lsq/__init__.py

@@ -0,0 +1,7 @@
+"""Module contains last-squares algorithms."""


Typo: last-squares. Also "This module".

ev-br · 2015-09-16T22:46:17Z

About huge sparse matrices: I just couldn't test in on my machine, as I said, it is unlikely to create troubles > (time / our of ram issues will come first).

Yeah, it's quite hard to believe it'll cause much trouble. FWIW, I'd punt on this.

So, it seems the only thing which holds this PR is a docstring example :-)

ev-br · 2015-09-21T11:39:50Z

OK, I keep maintaining that the last docstring example should be moved to the tutorial, but this can be done in a follow-up PR for the item 3 of #5044 (comment).

ev-br · 2015-09-21T11:48:25Z

Time to merge this, I'd think. (Further review is always welcome of course).

Thank you Nikolay, it's a great one.

[GSoC] ENH: New least-squares algorithms

nmayorov · 2015-09-21T12:11:35Z

Sorry Evgeni, that I was so slow for making final changes. But great to see it merged!

charris · 2015-09-21T12:53:15Z

Yay!

ev-br · 2015-10-17T15:48:20Z

Now that curve-fit PR is in, the plan #5044 (comment) is done save item 3 "docs/tutorial" :-).

nmayorov · 2015-10-25T19:50:15Z

Hi, Evgeni! Can you explain to me what is the current status of using IPython Notebooks as examples in scipy documentation?

ev-br · 2015-10-26T11:00:38Z

Nice to see you back Nikolay :-). The discussion is going on in #5233: AFAIU, the only last thing to agree upon is whether to host them in a separate repository or add them to the main scipy repo.

I think the best way to get the ball rolling is to send a WIP PR with your tweaks to the tutorial and links to your notebooks with examples.

anntzer · 2016-01-09T08:57:09Z

scipy/optimize/_lsq/least_squares.py

+    fun : callable
+        Function which computes the vector of residuals with the signature
+        ``fun(x, *args, **kwargs)``, i.e., the minimization proceeds with
+        respect to it's first argument. The argument ``x`` passed to this


typo: "it's" -> "its"

Thanks Antony --- would you care to send a PR with the fix please
09.01.2016 11:57 пользователь "Antony Lee" notifications@github.com
написал:

In scipy/optimize/_lsq/least_squares.py
#5044 (comment):

We call f(x) as a vector of residuals or simply residuals, and F(x) as a

cost function or simply cost.

We call rho(s) as a loss function, its purpose to reduce the influence

of outliers on the solution.

Partial derivatives of f with respect to x form m-by-n matrix called

Jacobian, where an element (i, j) equals the partial derivative of f[i]

with respect to x[j].

Parameters

fun : callable

Function which computes the vector of residuals with the signature

`fun(x, *args, **kwargs)`, i.e., the minimization proceeds with

respect to it's first argument. The argument `x` passed to this

typo: "it's" -> "its"

—
Reply to this email directly or view it on GitHub
https://github.com/scipy/scipy/pull/5044/files#r49262246.

Any corrections to the text would be helpful. Likely there are more mistakes and generally bad usage of English.

ev-br added enhancement A new feature or improvement scipy.optimize labels Jul 13, 2015

ev-br mentioned this pull request Jul 13, 2015

least_squares interface #5020

Closed

nmayorov force-pushed the bounded_lsq_sparse branch 2 times, most recently from 4ca0edb to 7f4e942 Compare July 19, 2015 01:18

ev-br reviewed Jul 22, 2015
View reviewed changes

nmayorov force-pushed the bounded_lsq_sparse branch from bc58791 to 36a0db4 Compare July 22, 2015 18:05

ev-br reviewed Jul 22, 2015
View reviewed changes

nmayorov force-pushed the bounded_lsq_sparse branch from 36a0db4 to 789075d Compare July 22, 2015 18:11

ev-br reviewed Jul 22, 2015
View reviewed changes

ev-br mentioned this pull request Jul 22, 2015

[GSoC] WIP: Bounded LSQ algorithms #5019

Closed

ev-br reviewed Jul 22, 2015
View reviewed changes

nmayorov force-pushed the bounded_lsq_sparse branch from 3d9370e to 7aabb06 Compare July 24, 2015 07:07

charris reviewed Sep 9, 2015
View reviewed changes

STY: Several stylistic improvements in optimize._lsq

b3766ad

ev-br mentioned this pull request Sep 14, 2015

wishlist: root finding with box constraints #5249

Closed

larsmans reviewed Sep 15, 2015
View reviewed changes

DOC: Fixed some errors/typos in docstrings of optimize._lsq

104cc8c

ev-br added this to the 0.17.0 milestone Sep 16, 2015

ev-br added a commit that referenced this pull request Sep 21, 2015

Merge pull request #5044 from nmayorov/bounded_lsq_sparse

1982080

[GSoC] ENH: New least-squares algorithms

ev-br merged commit 1982080 into scipy:master Sep 21, 2015

ev-br mentioned this pull request Sep 21, 2015

No way to set ranges for fitting parameters in Optimize functions #3129

Closed

ev-br mentioned this pull request Sep 21, 2015

Class-Based Levenberg-Marquardt Fitter #90

Closed

nmayorov deleted the bounded_lsq_sparse branch October 4, 2015 01:35

newville mentioned this pull request Nov 23, 2015

Fitting fails when 'min=' attribute is included for one parameter lmfit/lmfit-py#278

Closed

anntzer reviewed Jan 9, 2016
View reviewed changes

		return groups


		def group_sparse(int m, int n, int [:] indices, int [:] indptr):

		@@ -0,0 +1,7 @@
		"""Module contains last-squares algorithms."""

[GSoC] ENH: New least-squares algorithms #5044

[GSoC] ENH: New least-squares algorithms #5044

Conversation

nmayorov commented Jul 13, 2015

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ev-br commented Jul 22, 2015

nmayorov commented Jul 23, 2015

ev-br commented Jul 23, 2015

ev-br commented Jul 23, 2015

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

nmayorov commented Sep 13, 2015

ev-br commented Sep 14, 2015

nmayorov commented Sep 15, 2015

Choose a reason for hiding this comment

ev-br commented Sep 16, 2015

ev-br commented Sep 21, 2015

ev-br commented Sep 21, 2015

nmayorov commented Sep 21, 2015

charris commented Sep 21, 2015

ev-br commented Oct 17, 2015

nmayorov commented Oct 25, 2015

ev-br commented Oct 26, 2015

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment