Numpy/Scipyで二重対角行列や三重対角行列を作成

投稿日: 2019年1月31日 作成者: Funmatu

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
)

e38380e382a6e383b3e383ade383bce383895

download001

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

e38380e382a6e383b3e383ade383bce383894

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_diags

この方向性で考えるのはあまり良くないかも．

前述のように，「extend_eye」で考えるのが良さそう．

　
　

関連：
行列の対角に任意のベクトルを穴埋めする様な処理 – 対角線上に任意の配列、それ以外の場所に0を含む2次元配列を作成

numpy.ndarrayの各対角線の合計を求める

　

　

試行錯誤メモ（ゴミ箱の中身的なもの）．

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
)

try001

カテゴリー: 未分類 パーマリンク

Numpy/Scipyで二重対角行列や三重対角行列を作成 への2件のフィードバック

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

  2. ピンバック: 対角行列 | 粉末@それは風のように (日記)

コメントは受け付けていません。