numpy.array + scipy.sparse = numpy.matrixになる件

大規模データを扱うときnumpy.arrayのみではメモリに乗り切らないことがあります。
そんなとき役立つのがscipy.sparse(csr_matrix, csc_matrix)なのですが、それらを用いた演算をしたとき何故かそれまで使っていないnumpy.matrixなるデータ型が返ってきました。   
この記事はその事象について調べた結果を自分用メモとしてまとめたものです（正直かなりの部分[scipyのissue](https://github.com/scipy/scipy/issues/7510)に載っている内容なのですが）。

### TL;DR
- array +-/* sparse = matrixという挙動は正常
- 正直戻り値をmatrixでない別物にしたいようだが修正にはかなり時間が掛かる模様
- sparseとの演算結果をsparseのみで行うにはすごく工夫が必要そう。

### 背景
協調フィルタリングの文脈で評価値行列など諸々をndarrayで持っていたのですが、とてもではないけれどもメモリに収まりきらずsparse matrixを用いることになりました。
sparse matrixを用いて類似度を計算していたところ、以下のような現象に遭遇しました。

```python
# sp_matrix ... sparse.csr_matrix
# np_matrix ... numpy.array

test_ matrix = np_matrix * sp_matrix
-> array * sparse_matrixのはずなのに戻り値がnumpy.matrixになる！
```

そもそもnumpy.matrixなんて使ったことなかったのでそもそもそれってなんなの？となり、しかもなんで戻り値が変わるねんと疑問に感じました。
以上が事の顛末です。

### そもそも
この辺は知ってる方からしたらおさらいです。

#### 1. numpy.array
みなさんおなじみのndarrayです。
通常機械学習とかで行列扱うときはこれを使うのではないでしょうか。

[公式ドキュメント](https://docs.scipy.org/doc/numpy/reference/generated/numpy.array.html)より
> An array, any object exposing the array interface, an object whose __array__ method returns an array, or any (nested) sequence.


#### 2. numpy.matrix
一方聞き馴染みがないのはこちらのnumpy.matrixです。
こちらも[公式ドキュメント](https://docs.scipy.org/doc/numpy-1.15.1/reference/generated/numpy.matrix.html)からどういうものか確認してみましょう。
> Returns a matrix from an array-like object, or from a string of data. A matrix is a specialized 2-D array that retains its 2-D nature through operations. It has certain special operators, such as * (matrix multiplication) and ** (matrix power).

なにやら基本的に2次元配列を扱うものらしいです。
しかしこの公式ドキュメント上部にはこのような記載があります。

> Note
> It is no longer recommended to use this class, even for linear algebra. Instead use regular arrays. The class may be removed in the future.

要するに今は使わないでくれってことですね。


#### 3. scipy.sparse.xxx_matrix

こちらも使っている方にはおなじみのsparse matrixです。疎行列を扱うときに使われます。
xxx_matrixとしているのは色々種類があるからで、用途に応じてlil_matrix, csc_matrix, csr_matrixといったものが用意されています。

例によって以下[公式ドキュメント](https://docs.scipy.org/doc/scipy/reference/sparse.html)ですが、ちょっと端的すぎてわからんです。
> SciPy 2-D sparse matrix package for numeric data.

使い分け等々に関して私は[こちら](https://stackoverflow.com/questions/36969886/using-a-sparse-matrix-versus-numpy-array)の記事を参考にしました（のですみませんがこの辺の説明は端折ります）。


### 実験

それぞれがどういうものかわかったところでここから実験です。

今回は以下の方々を用いて実験します。
csr_matrixでやっていますが、元のissueでもやられているようにcsc_matrixでも同様の結果になります。

```python
import numpy as np
from scipy.sparse import csr_matrix, csc_matrix

data = [[1,1,0],
        [1,0,0],
        [0,1,0]]

>>> array = np.array(data)
>>> print(array)
[[1 1 0]
 [1 0 0]
 [0 1 0]]
>>> matrix = np.matrix(data)
>>> print(matrix)
[[1 1 0]
 [1 0 0]
 [0 1 0]]
>>> sparse = csr_matrix(data)
>>> print(sparse)
  (0, 0)	1
  (0, 1)	1
  (1, 0)	1
  (2, 1)	1
```

ここからそれぞれのケースを実験していきます。
ちなみにここでやっていることは[こちらの参考リンク](https://github.com/scipy/scipy/issues/7510)を概ね模倣したものです。

#### 1. array + sparse
```python
>>> array + sparse
matrix([[2, 2, 0],
        [2, 0, 0],
        [0, 2, 0]])

>>> array - sparse
matrix([[0, 0, 0],
        [0, 0, 0],
        [0, 0, 0]], dtype=int64)

>>> array * sparse
array([[2, 1, 0],
       [1, 1, 0],
       [1, 0, 0]], dtype=int64)
```
足し引きの方はmatrixで返ってきており、積の方だけarrayで返ってきています。

#### 2. matrix + sparse

```python
>>> matrix + sparse
matrix([[2, 2, 0],
        [2, 0, 0],
        [0, 2, 0]])
>>> matrix - sparse
matrix([[0, 0, 0],
        [0, 0, 0],
        [0, 0, 0]], dtype=int64)
>>> matrix * sparse
matrix([[2, 1, 0],
        [1, 1, 0],
        [1, 0, 0]], dtype=int64)
```
こちらの場合戻り値は全てmatrixです。

### 結論
#### 1. 基本的にnumpy.array +-\* scipy.sparse = numpy.matrixが意図した挙動。
  これはscipy.sparseのAPI設計がnumpy.matrixを模したものであることがその由来。
  問題として捉えてはいるものの修正には時間が掛かるようです。
  念の為原文（一部）を転記します。


> In general, I agree that users should avoid numpy.matrix whenever possible.
> Unfortunately, the scipy.sparse API design decision to mimic numpy.matrix can't be changed now.
> There's a growing list of efforts to produce a sparse ndarray-like library that may someday make it into scipy, but that will take considerable time.


#### 2. array * sparse = array, およびarray += sparseはバグ
  ただし上述のような事情のため修正には時間がかかりそうです。


### 雑感
個人的にはそもそも{array, matrix} * sparseの結果がsparseであってもdenseな結果が返ってくるというのはちょっとした発見でした。
類似度計算の途中どうしてもarrayができてしまう[^1]ので、「array * sparseで結果もsparse」になってくれないかなーと思いましたがそうは行かないようです。
リンク先のissueでも言われているように、そのarray * sparseの結果がsparseな行列であったとしても型としてはnumpy.arrayになってしまう…バグでもそれとは辛い。

ということで以上です。

### 参考リンク
- [Operations between numpy.array and scipy.sparse matrix return inconsistent array type #7510](https://github.com/scipy/scipy/issues/7510)   
←今回の記事に最も関係するissue   
- [Compute all pairwise vector similarities within a sparse matrix (Python)](http://na-o-ys.github.io/others/2015-11-07-sparse-vector-similarities.html)
- [Compute Jaccard distances on sparse matrix](https://stackoverflow.com/questions/32805916/compute-jaccard-distances-on-sparse-matrix)   
←こちら2つはjaccard指数をsparse matrixで出すとき参考にしたブログとやり取り

[^1]: 例えば(sparse).getnnz()や(sparse).diagonal()みたいなやつも戻り値はarrayになる。