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

大規模データを扱うときnumpy.arrayのみではメモリに乗り切らないことがあります。
そんなとき役立つのがscipy.sparse(csr_matrix, csc_matrix)なのですが、それらを用いた演算をしたとき何故かそれまで使っていないnumpy.matrixなるデータ型が返ってきました。

この記事はその事象について調べた結果を自分用メモとしてまとめたものです（正直かなりの部分scipyのissueに載っている内容なのですが）。

TL;DR

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

背景

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

# 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です。
通常機械学習とかで行列扱うときはこれを使うのではないでしょうか。

公式ドキュメントより

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です。
こちらも公式ドキュメントからどういうものか確認してみましょう。

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といったものが用意されています。

例によって以下公式ドキュメントですが、ちょっと端的すぎてわからんです。

SciPy 2-D sparse matrix package for numeric data.

使い分け等々に関して私はこちらの記事を参考にしました（のですみませんがこの辺の説明は端折ります）。

実験

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

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

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

ここからそれぞれのケースを実験していきます。
ちなみにここでやっていることはこちらの参考リンクを概ね模倣したものです。

1. array + sparse

>>> 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

>>> 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ができてしまう¹ので、「array * sparseで結果もsparse」になってくれないかなーと思いましたがそうは行かないようです。
リンク先のissueでも言われているように、そのarray * sparseの結果がsparseな行列であったとしても型としてはnumpy.arrayになってしまう…バグでもそれとは辛い。

ということで以上です。

参考リンク

Operations between numpy.array and scipy.sparse matrix return inconsistent array type #7510
←今回の記事に最も関係するissue
Compute all pairwise vector similarities within a sparse matrix (Python)
Compute Jaccard distances on sparse matrix
←こちら2つはjaccard指数をsparse matrixで出すとき参考にしたブログとやり取り

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