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を用いて行列内にベクトルを埋め込む – 行列の対角に任意の行列を穴埋めする様な処理 | 粉末@それは風のように (日記)
ピンバック: 対角行列 | 粉末@それは風のように (日記)