Repository: fatiando/verde

Issue:
Add optional argument to avoid overwriting Jacobian in least square
**Description of the desired feature**

The `verde.base.least_squares` solves a linear problem with Tikhonov zero order regularization by rescaling the Jacobian matrix to unit variance columns, so the damping parameter is independent of the units of the misfit. While scaling the Jacobian, a `copy=False` argument is passed to `StandardScaler` so we can save memory by scaling the Jacobian in place. This has the advantage of saving memory, which is a good feature to have for most cases. But on some examples we may want to reuse the computed Jacobian, so we want to avoid to overwrite it.

Would be nice to have an optional `copy` argument on `verde.base.least_squares` that is passed to the `StandardScaler` constructor.

If we decide not  to add this option, I would suggest to add a warning on the docstring explicitly saying that the Jacobian matrix will be overwritten. Also to change the following line:
https://github.com/fatiando/verde/blob/9884a274fc7cc49e213345d21f17d960082eaef0/verde/base/least_squares.py#L48

to:
```python
scaler.fit_transform(jacobian)
```
because `fit_transform` scales the matrix it in place if `copy` is `False`.
This way we can also figure out that the matrix will be overwritten by just reading the code.

**Are you willing to help implement and maintain this feature?** Yes


Discussion:

Patch:
diff --git a/verde/base/least_squares.py b/verde/base/least_squares.py
index 81ec3b3..e2b738a 100644
--- a/verde/base/least_squares.py
+++ b/verde/base/least_squares.py
@@ -7,7 +7,7 @@ from sklearn.preprocessing import StandardScaler
 from sklearn.linear_model import LinearRegression, Ridge
 
 
-def least_squares(jacobian, data, weights, damping=None):
+def least_squares(jacobian, data, weights, damping=None, copy_jacobian=False):
     """
     Solve a weighted least-squares problem with optional damping regularization
 
@@ -17,6 +17,12 @@ def least_squares(jacobian, data, weights, damping=None):
     required for predictions. Doesn't normalize the column means because that
     operation can't be undone.
 
+    .. warning::
+
+        Setting `copy_jacobian` to True will copy the Jacobian matrix, doubling
+        the memory required. Use it only if the Jacobian matrix is needed
+        afterwards.
+
     Parameters
     ----------
     jacobian : 2d-array
@@ -31,6 +37,9 @@ def least_squares(jacobian, data, weights, damping=None):
     damping : None or float
         The positive damping (Tikhonov 0th order) regularization parameter. If
         ``damping=None``, will use a regular least-squares fit.
+    copy_jacobian: bool
+        If False, the Jacobian matrix will be scaled inplace. If True, the
+        Jacobian matrix will be copied before scaling. Default False.
 
     Returns
     -------
@@ -44,7 +53,7 @@ def least_squares(jacobian, data, weights, damping=None):
                 jacobian.shape
             )
         )
-    scaler = StandardScaler(copy=False, with_mean=False, with_std=True)
+    scaler = StandardScaler(copy=copy_jacobian, with_mean=False, with_std=True)
     jacobian = scaler.fit_transform(jacobian)
     if damping is None:
         regr = LinearRegression(fit_intercept=False, normalize=False)


Test patch:
diff --git a/verde/tests/test_base.py b/verde/tests/test_base.py
index 4f2f61d..88ba369 100644
--- a/verde/tests/test_base.py
+++ b/verde/tests/test_base.py
@@ -6,6 +6,7 @@ import numpy as np
 import numpy.testing as npt
 import pytest
 
+from ..base.least_squares import least_squares
 from ..base.utils import check_fit_input, check_coordinates
 from ..base.base_classes import (
     BaseGridder,
@@ -232,3 +233,16 @@ def test_check_fit_input_fails_weights():
         check_fit_input(coords, data, weights)
     with pytest.raises(ValueError):
         check_fit_input(coords, (data, data), weights)
+
+
+def test_least_squares_copy_jacobian():
+    """
+    Test if Jacobian matrix is copied or scaled inplace
+    """
+    jacobian = np.identity(5)
+    original_jacobian = jacobian.copy()
+    data = np.array([1, 2, 3, 4, 5], dtype=float)
+    least_squares(jacobian, data, weights=None, copy_jacobian=True)
+    npt.assert_allclose(jacobian, original_jacobian)
+    least_squares(jacobian, data, weights=None)
+    assert not np.allclose(jacobian, original_jacobian)