Rolling window viewを用いた効率的な計算

Efficient way to lag vector multiplication using numpy – StackOverflow

まず,OPのを愚直に求めてみると,

import numpy as np


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


a = np.array([3, 2, 0])
b = np.array([2, 5, 6, 3, 8, 4])
a_ = extend_eye(a, 6)[:4]
b_ = np.diag(b)
res = (a_[:, None]*b_).sum(-1)

display(a_)
display(b_)
res
array([[3, 2, 0, 0, 0, 0],
       [0, 3, 2, 0, 0, 0],
       [0, 0, 3, 2, 0, 0],
       [0, 0, 0, 3, 2, 0]])
array([[2, 0, 0, 0, 0, 0],
       [0, 5, 0, 0, 0, 0],
       [0, 0, 6, 0, 0, 0],
       [0, 0, 0, 3, 0, 0],
       [0, 0, 0, 0, 8, 0],
       [0, 0, 0, 0, 0, 4]])
array([[ 6, 10,  0,  0,  0,  0],
       [ 0, 15, 12,  0,  0,  0],
       [ 0,  0, 18,  6,  0,  0],
       [ 0,  0,  0,  9, 16,  0]])

形状に意味がある訳ではなく,値が欲しい様なので,そうするとベクトルaの形に合わせて,ベクトルbのrolling windowを求めて乗算すれば良い.

def f(a, N):
    s = a.strides[0]
    strided = np.lib.stride_tricks.as_strided
    return strided(a, shape=(len(a)-N+1, N), strides=(s, s))


f(b, len(a)) * a
array([[ 6, 10,  0],
       [15, 12,  0],
       [18,  6,  0],
       [ 9, 16,  0]])
カテゴリー: 未分類 パーマリンク

コメントを残す

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

WordPress.com ロゴ

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

Google フォト

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

Twitter 画像

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

Facebook の写真

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

%s と連携中

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