Create a bidiagonal matrix – StackOverflow
OPで提示される様な単純なケースでは,OPが示している方法は決して悪くないと思う.
import numpy as np main_diag = [1,2,3,4,5] off1 = [1,2,3,4] np.diag(main_diag) + np.diag(off1, 1)
array([[1, 1, 0, 0, 0], [0, 2, 2, 0, 0], [0, 0, 3, 3, 0], [0, 0, 0, 4, 4], [0, 0, 0, 0, 5]])
ただ,これはハードコーディングなので,汎用性がない(必ずしもそれが悪い訳では無いと思うけど).なので,汎用性がある方法は?というのが質問の趣旨だろう.この手の処理をしたい場合,numpy.lib.stride_tricks.as_stridedを用いるのがリーズナブルチョイス.任意のN重対角行列を求めたい場合,関連のextend_eye関数(引数はnumpy.eyeと対応付けているつもり)の様な方法を用いれば良い.
def extend_eye(p, N, M=None, k=0, dtype=None): if M is None: M = N if dtype is None: dtype = p.dtype n = N - 1 b = np.zeros((n+M,), dtype=dtype) b[n+k:n+len(p)+k] = p s = b.strides[0] strided = np.lib.stride_tricks.as_strided return strided(b[n:], shape=(N, M), strides=(-s, s)) N = 2 p = np.array([1]*N) extend_eye(p, 5) * np.arange(1, 6)[:, None]
array([[1, 1, 0, 0, 0], [0, 2, 2, 0, 0], [0, 0, 3, 3, 0], [0, 0, 0, 4, 4], [0, 0, 0, 0, 5]])
関連のextend_eye関数を流用しているので,この解法に最適化されていないが,汎用性は高いし,回答にあるscipy.sparse.diagsと比較しても,(対角数が任意の場合)最大で20倍位速い.2重対角行列の様なケースでも,N数が増加してもほぼ同じである事を考えれば,こちらの方が十分に効率的だしリーズナブル(そして,結果的にあれこれ悩むよりもそれが一番良かった).
import numpy as np from scipy import sparse import perfplot def extend_eye(p, N, M=None, k=0, dtype=None): if M is None: M = N if dtype is None: dtype = p.dtype n = N - 1 b = np.zeros((n+M,), dtype=dtype) b[n+k:n+len(p)+k] = p s = b.strides[0] strided = np.lib.stride_tricks.as_strided return strided(b[n:], shape=(N, M), strides=(-s, s)) def sparse_diags(n, N): diags = [range(1, N+1-i) for i in range(n)] return sparse.diags(diags, range(n)).toarray() def myfunc(n, N): p = np.array([1]*n) return extend_eye(p, N) * np.arange(1, N+1)[:, None] perfplot.show( setup=lambda N: [2, N], kernels=[ lambda x: myfunc(x[0], x[1]), lambda x: sparse_diags(x[0], x[1]) ], labels=['extend_eye', 'scipy.sparse.diags'], n_range=[2**k for k in range(1, 16)], xlabel='N', logx=True, logy=True, equality_check=False )
%timeit sparse_diags(2, 10) %timeit myfunc(2, 10) %timeit sparse_diags(2, 100) %timeit myfunc(2, 100) %timeit sparse_diags(2, 1000) %timeit myfunc(2, 1000) %timeit sparse_diags(2, 10000) %timeit myfunc(2, 10000)
1000 loops, best of 3: 244 µs per loop 100000 loops, best of 3: 17.3 µs per loop 1000 loops, best of 3: 291 µs per loop 10000 loops, best of 3: 40.5 µs per loop 10 loops, best of 3: 35.6 ms per loop 10 loops, best of 3: 35.5 ms per loop 1 loop, best of 3: 242 ms per loop 1 loop, best of 3: 253 ms per loop
import numpy as np from scipy import sparse import perfplot import random def extend_eye(p, N, M=None, k=0, dtype=None): if M is None: M = N if dtype is None: dtype = p.dtype n = N - 1 b = np.zeros((n+M,), dtype=dtype) b[n+k:n+len(p)+k] = p s = b.strides[0] strided = np.lib.stride_tricks.as_strided return strided(b[n:], shape=(N, M), strides=(-s, s)) def sparse_diags(n, N): diags = [range(1, N+1-i) for i in range(n)] return sparse.diags(diags, range(n)).toarray() def myfunc(n, N): p = np.array([1]*n) return extend_eye(p, N) * np.arange(1, N+1)[:, None] perfplot.show( setup=lambda N: [random.randrange(1, N), N], kernels=[ lambda x: myfunc(x[0], x[1]), lambda x: sparse_diags(x[0], x[1]) ], labels=['extend_eye', 'scipy.sparse.diags'], n_range=[2**k for k in range(1, 16)], xlabel='N', logx=True, logy=True, equality_check=False )
scipy.sparse.diagsは大前提として(当然ながら)目的のアレイがスパースである場合でかつ,オフセットが乱雑なケースを想定する場合,簡便かも知れないが,DIAgonal formatで効率的なフォーマットではないので,CSRやCSCに置き換えない限りは空間効率も特に良くはないし,時間効率は良くないので,必ずしもリーズナブルだとは思わない(より何かしらの仮定が置ける場合,スパース行列として考える合理性).
ちなみに,コメント欄でフラットフラット言っているのは,多分こういう事だろう.
import numpy as np main_diag = [1,2,3,4,5] a = np.diag(main_diag) idx = np.flatnonzero(a) a.flat[idx[:-1]+1] = np.arange(1, 5) a
array([[1, 1, 0, 0, 0], [0, 2, 2, 0, 0], [0, 0, 3, 3, 0], [0, 0, 0, 4, 4], [0, 0, 0, 0, 5]])
ただ,これを一般化しようと思うと,はっきり言って面倒くさいと思う(というか,これを一般化すると,「extend_eye」の様な処理になるのでは).
numpy.lib.stride_tricks.as_stridedを用いた解法を,このケースにフィッティングしようと思うと,パターンベクトルをブロードキャストして欲しい形状にした後,両端にゼロベクトルを結合して,ローリングされたビューを求めれば良さそう.……と思ったけど,巧いこと思い付かなかった.例えば,
import numpy as np def extend_diags(n, N, dtype=int): p = np.broadcast_to(np.arange(1, N+1, dtype=dtype)[:, None], (N, n)) a = np.concatenate((np.zeros((N, N), dtype=dtype), p, np.zeros((N, N-n+1), dtype=dtype)), 1) s1, s2 = a.strides strided = np.lib.stride_tricks.as_strided return strided(a, shape=(N+1, N*3), strides=(s1+s2, -s2))[1:, N-n+3:2*N-n+3] print(extend_diags(2, 5)) print(extend_diags(2, 10)) print(extend_diags(3, 5)) print(extend_diags(3, 10)) print(extend_diags(4, 5)) print(extend_diags(4, 10))
[[1 1 0 0 0] [0 2 2 0 0] [0 0 3 3 0] [0 0 0 4 4] [0 0 0 0 5]] [[ 1 1 0 0 0 0 0 0 0 0] [ 0 2 2 0 0 0 0 0 0 0] [ 0 0 3 3 0 0 0 0 0 0] [ 0 0 0 4 4 0 0 0 0 0] [ 0 0 0 0 5 5 0 0 0 0] [ 0 0 0 0 0 6 6 0 0 0] [ 0 0 0 0 0 0 7 7 0 0] [ 0 0 0 0 0 0 0 8 8 0] [ 0 0 0 0 0 0 0 0 9 9] [ 0 0 0 0 0 0 0 0 0 10]] [[1 1 1 0 0] [0 2 2 2 0] [0 0 3 3 3] [0 0 0 4 4] [0 0 0 0 5]] [[ 1 1 1 0 0 0 0 0 0 0] [ 0 2 2 2 0 0 0 0 0 0] [ 0 0 3 3 3 0 0 0 0 0] [ 0 0 0 4 4 4 0 0 0 0] [ 0 0 0 0 5 5 5 0 0 0] [ 0 0 0 0 0 6 6 6 0 0] [ 0 0 0 0 0 0 7 7 7 0] [ 0 0 0 0 0 0 0 8 8 8] [ 0 0 0 0 0 0 0 0 9 9] [ 0 0 0 0 0 0 0 0 0 10]] [[1 1 1 1 0] [0 2 2 2 2] [0 0 3 3 3] [0 0 0 4 4] [0 0 0 0 5]] [[ 1 1 1 1 0 0 0 0 0 0] [ 0 2 2 2 2 0 0 0 0 0] [ 0 0 3 3 3 3 0 0 0 0] [ 0 0 0 4 4 4 4 0 0 0] [ 0 0 0 0 5 5 5 5 0 0] [ 0 0 0 0 0 6 6 6 6 0] [ 0 0 0 0 0 0 7 7 7 7] [ 0 0 0 0 0 0 0 8 8 8] [ 0 0 0 0 0 0 0 0 9 9] [ 0 0 0 0 0 0 0 0 0 10]]
ただ,これは時空間効率が(明らかに)悪い.
import numpy as np from scipy import sparse import perfplot def extend_eye(p, N, M=None, k=0, dtype=None): if M is None: M = N if dtype is None: dtype = p.dtype n = N - 1 b = np.zeros((n+M,), dtype=dtype) b[n+k:n+len(p)+k] = p s = b.strides[0] strided = np.lib.stride_tricks.as_strided return strided(b[n:], shape=(N, M), strides=(-s, s)) def sparse_diags(n, N): diags = [range(1, N+1-i) for i in range(n)] return sparse.diags(diags, range(n)).toarray() def myfunc(n, N): p = np.array([1]*n) return extend_eye(p, N) * np.arange(1, N+1)[:, None] def extend_diags(n, N, dtype=int): p = np.broadcast_to(np.arange(1, N+1, dtype=dtype)[:, None], (N, n)) a = np.concatenate((np.zeros((N, N), dtype=dtype), p, np.zeros((N, N-n+1), dtype=dtype)), 1) s1, s2 = a.strides strided = np.lib.stride_tricks.as_strided return strided(a, shape=(N+1, N*3), strides=(s1+s2, -s2))[1:, N-n+3:2*N-n+3] def two_diags(N, dtype=int): n = 2 N += 1 p = np.broadcast_to(np.arange(N, dtype=dtype)[:, None], (N, n)) a = np.concatenate((p, np.zeros((N, N), dtype=dtype)), 1) s1, s2 = a.strides strided = np.lib.stride_tricks.as_strided return strided(a, shape=(N, N+n), strides=(s1+s2, -s2))[1:, :N-1] perfplot.show( setup=lambda N: [2, N], kernels=[ lambda x: myfunc(x[0], x[1]), lambda x: sparse_diags(x[0], x[1]), lambda x: extend_diags(x[0], x[1]), lambda x: two_diags(x[1]), ], labels=['extend_eye x broadcast', 'scipy.sparse.diags', 'extend_diags', 'two_diags'], n_range=[2**k for k in range(1, 15)], xlabel='N', logx=True, logy=True, equality_check=False )
この方向性で考えるのはあまり良くないかも.
前述のように,「extend_eye」で考えるのが良さそう.
関連:
行列の対角に任意のベクトルを穴埋めする様な処理 – 対角線上に任意の配列、それ以外の場所に0を含む2次元配列を作成
試行錯誤メモ(ゴミ箱の中身的なもの).
import numpy as np from scipy import sparse import perfplot def extend_eye(p, N, M=None, k=0, dtype=None): if M is None: M = N if dtype is None: dtype = p.dtype n = N - 1 b = np.zeros((n+M,), dtype=dtype) b[n+k:n+len(p)+k] = p s = b.strides[0] strided = np.lib.stride_tricks.as_strided return strided(b[n:], shape=(N, M), strides=(-s, s)) def sparse_diags(n, N): diags = [range(1, N+1-i) for i in range(n)] return sparse.diags(diags, range(n)).toarray() def myfunc(n, N): p = np.array([1]*n) return extend_eye(p, N) * np.arange(1, N+1)[:, None] def extend_diags(n, N, dtype=np.uint16): res = np.ones((N, N), dtype=dtype) res *= np.arange(1, N+1, dtype=dtype)[:, None] res[np.mask_indices(N, np.triu, n)] = 0 res[np.mask_indices(N, np.tril, -1)] = 0 return res def two_diags(N, dtype=int): n = 2 res = np.zeros((N, N), dtype=dtype) ridx, cidx = np.diag_indices(N) idx = ((ridx[:-1]*N + cidx[:-1])[:, None] + np.arange(n)).ravel() res.flat[idx] = (ridx+1).repeat(n) res.flat[-1] = N return res perfplot.show( setup=lambda N: [2, N], kernels=[ lambda x: myfunc(x[0], x[1]), lambda x: sparse_diags(x[0], x[1]), lambda x: extend_diags(x[0], x[1]), lambda x: two_diags(x[1]) ], labels=['extend_eye x broadcast', 'scipy.sparse.diags', 'numpy_mask_indices', 'flat indexing'], n_range=[2**k for k in range(1, 15)], xlabel='N', logx=True, logy=True, equality_check=False )
ピンバック: slicing/indexingを用いて行列内にベクトルを埋め込む – 行列の対角に任意の行列を穴埋めする様な処理 | 粉末@それは風のように (日記)
ピンバック: 対角行列 | 粉末@それは風のように (日記)