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

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. ピンバック: 対角行列 | 粉末@それは風のように (日記)

コメントを残す

以下に詳細を記入するか、アイコンをクリックしてログインしてください。

WordPress.com ロゴ

WordPress.com アカウントを使ってコメントしています。 ログアウト /  変更 )

Google フォト

Google アカウントを使ってコメントしています。 ログアウト /  変更 )

Twitter 画像

Twitter アカウントを使ってコメントしています。 ログアウト /  変更 )

Facebook の写真

Facebook アカウントを使ってコメントしています。 ログアウト /  変更 )

%s と連携中

このサイトはスパムを低減するために Akismet を使っています。コメントデータの処理方法の詳細はこちらをご覧ください