Repository: sympy/sympy
Pull request: give SparseMatrix own reshape routine
<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234" (see
https://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->
closes #13828

#### Brief description of what is fixed or changed


#### Other comments


#### Release Notes

<!-- Write the release notes for this release below between the BEGIN and END
statements. The basic format is a bulleted list with the name of the subpackage
and the release note for this PR. For example:

* solvers
  * Added a new solver for logarithmic equations.

* functions
  * Fixed a bug with log of integers. Formerly, `log(-x)` incorrectly gave `-log(x)`.

* physics.units
  * Corrected a semantical error in the conversion between volt and statvolt which
    reported the volt as being larger than the statvolt.

or if no release note(s) should be included use:

NO ENTRY

See https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more
information on how to write release notes. The bot will check your release
notes automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
* matrices
  * SparseMatrix reshaping is now more efficient
<!-- END RELEASE NOTES -->


Problem statement:
SparseMatrix should have its own reshape routine
The following should be about a million times faster with a dedicated reshape routine that just moves around existing elements and does nothing with 0s:

```
>>> s = SparseMatrix.zeros(1000, 1000)
>>> s[0] = 1
>>> s.reshape(1, 10**6)
```

Hints from issue discussion:
@smichr I want to work on this issue. I need some help that routine should create another zero matrix and substitute non zero element or it should modify existing one?
The same thing happens with this
`>>> s = eye(10**5)`
`Traceback (most recent call last):`
`
File "<string>", line 1, in <module>`
`  File "/base/data/home/apps/s~sympy-live-hrd/49.400254913747479351/sympy/sympy/matrices/dense.py", line 1275, in eye`
    `return cls.eye(n)`
  `File "/base/data/home/apps/s~sympy-live-hrd/49.400254913747479351/sympy/sympy/matrices/dense.py", line 529, in eye`
  `  mat = [cls._sympify(0)]*n*n
MemoryError
`

Also, in this case, take very long time
`>>> M = eye(10**4)
`


should we try to improve some MemoryError so that it can be used for large data such as for machine learning problems and optimizations?

> The same thing happens with this
```
>>> s = eye(10**5)
Traceback (most recent call last):
 File "<string>", line 1, in <module>
  File "/base/data/home/apps/s~sympy-live-hrd/49.400254913747479351/sympy/sympy/matrices/dense.py", line 1275, in eye
return cls.eye(n)
File "/base/data/home/apps/s~sympy-live-hrd/49.400254913747479351/sympy/sympy/matrices/dense.py", line 529, in eye
 mat = [cls._sympify(0)]*n*n MemoryError
```
I wonder if we should have a FormulaMatrix (kind of like default dict) which uses a formula to return elements and doesn't autoexpand, so the following would take very little space yet logically correspond to `10**6` elements (like it does now) of the MutableDenseMatrix.
```
eye1000 = Matrix(10**3,10**3,lambda (i,j): 1 if i==j else 0)
```
There is already FunctionMatrix in the matrix expressions:

```py
>>> FunctionMatrix(10**3, 10**3, Lambda((i, j), Piecewise((1, Eq(i, j)), (0, True))))
FunctionMatrix(1000, 1000, Lambda((i, j), Piecewise((1, Eq(i, j)), (0, True))))
```

Of course, matrix expressions already have Identity directly:

```py
>>> Identity(10**3)
I
```

Matrix expressions are a much better approach to this problem because they can represent the expression symbolically, meaning it's possible to construct further matrices without evaluating. They are also unevaluated by default meaning you don't have to worry about an operation accidentally trying to evaluate a million matrix entries. 
The OP issue is not hard to fix:
```python
In [12]: M = zeros(1000, 1000)

In [13]: M[0,0] = 1

In [14]: M.todok()
Out[14]: {(0, 0): 1}

In [15]: Matrix.from_dok(*M.shape, M.todok()) == M
Out[15]: True
```
A reshape function can be added that uses `todok`. I haven't added it at the DomainMatrix level because reshaping a matrix isn't generally useful.

Commit message: give SparseMatrix own reshape routine

Patch to review:
diff --git a/sympy/matrices/matrixbase.py b/sympy/matrices/matrixbase.py
index 7080b19e8a03..3ac1c8c87697 100644
--- a/sympy/matrices/matrixbase.py
+++ b/sympy/matrices/matrixbase.py
@@ -505,7 +505,9 @@ def reshape(self, rows, cols):
         """
         if self.rows * self.cols != rows * cols:
             raise ValueError("Invalid reshape parameters %d %d" % (rows, cols))
-        return self._new(rows, cols, lambda i, j: self[i * cols + j])
+        dok = {divmod(i*self.cols + j, cols):
+            v for (i, j), v in self.todok().items()}
+        return self._eval_from_dok(rows, cols, dok)
 
     def row_del(self, row):
         """Delete the specified row."""
diff --git a/sympy/matrices/tests/test_sparse.py b/sympy/matrices/tests/test_sparse.py
index 4d257c8062f2..6b76c8caf07c 100644
--- a/sympy/matrices/tests/test_sparse.py
+++ b/sympy/matrices/tests/test_sparse.py
@@ -316,6 +316,7 @@ def sparse_zeros(n):
         SparseMatrix([(0, 1, 2), (3, 1, 2), (3, 4, 2), (3, 4, 5)])
     assert m1.reshape(2, 6) == \
         SparseMatrix([(0, 1, 2, 3, 1, 2), (3, 4, 2, 3, 4, 5)])
+    assert sparse_eye(10**3).reshape(1, 10**6).shape == (1, 10**6)  # is not slow
 
     # test_applyfunc
     m0 = sparse_eye(3)


Code review comments:
diff_hunk: @@ -316,6 +316,7 @@ def sparse_zeros(n):
         SparseMatrix([(0, 1, 2), (3, 1, 2), (3, 4, 2), (3, 4, 5)])
     assert m1.reshape(2, 6) == \
         SparseMatrix([(0, 1, 2, 3, 1, 2), (3, 4, 2, 3, 4, 5)])
+    assert sparse_eye(10**3).reshape(1, 10**6).shape == (1, 10**6)  # is not slow
path: sympy/matrices/tests/test_sparse.py
text: @user1:
I think it would be better just to remove this test. It isn't slow enough on master that it would get noticed in a PR. We can use benchmarks if we want to catch performance regressions.