연결된 행렬 곱셈

Question

조밀 한 행렬 A와 벡터 e, f에 대해 A * diag (e) * A^T * f의 빠른 곱셈에 대한 제안 사항은 무엇입니까?

이것은 내가 지금 가지고있는 것입니다.

v[:] = 0 
for i in range(N): 
    for j in range(N): 
     v[i] = v[i]+A[i,j]*e[j]*np.dot(A[:,j],f) 

감사합니다, @rubenvb에 의해

Answer 1

제안은 아마 그것을 할 수있는 가장 간단한 방법입니다. 또 다른 방법은 einsum을 사용하는 것입니다.

다음은 예입니다. 내가 사용하는거야 a, e 및 f 다음

In [95]: a 
Out[95]: 
array([[1, 2, 3], 
     [4, 5, 6], 
     [7, 8, 9]]) 

In [96]: e 
Out[96]: array([-1, 2, 3]) 

In [97]: f 
Out[97]: array([5, 4, 1]) 

이 NumPy와 코드로 수식의 직접 번역입니다.

: 당신은 그 인수와 관련된 인덱스 라벨을 교환하여 a 전치 할 필요성을 제거 할 수

In [99]: np.einsum('ij,j,jk,k', a, e, a.T, f) 
Out[99]: array([ 556, 1132, 1708]) 

:

여기

In [98]: a.dot(np.diag(e)).dot(a.T).dot(f) 
Out[98]: array([ 556, 1132, 1708]) 

가 einsum 버전 : 그것은 기본적으로 @의 rubenvb의 제안과 동일

In [100]: np.einsum('ij,j,kj,k', a, e, a, f) 
Out[100]: array([ 556, 1132, 1708]) 

Answer 2

@ rubenvb의 의견에 따르면 A.dot(np.diag(e)).dot(A.transpose()).dot(f)을 사용하는 것이 좋습니다. 그러나 실제로 2D 배열을 diag(e)으로 만들 필요가 없으므로 하나의 행렬 곱셈을 건너 뜁니다. 또한 A.T 및 f에 대한 장소를 바꿀 수 있으므로 전치도 방지 할 수 있습니다. 따라서, 단순하고 훨씬 더 효율적이 솔루션 같은 진화 것 -

A.dot(e*f.dot(A)) 

을 여기에 모든 제안 된 방법에 알맞은 크기의 배열에 빠른 런타임 테스트입니다 -

In [226]: # Setup inputs 
    ...: N = 200 
    ...: A = np.random.rand(N,N) 
    ...: e = np.random.rand(N,) 
    ...: f = np.random.rand(N,) 
    ...: 

In [227]: %timeit np.einsum('ij,j,kj,k', A, e, A, f) # @Warren Weckesser's soln 
10 loops, best of 3: 77.6 ms per loop 

In [228]: %timeit A.dot(np.diag(e)).dot(A.transpose()).dot(f) # @rubenvb's soln 
10 loops, best of 3: 18.6 ms per loop 

In [229]: %timeit A.dot(e*f.dot(A)) # Proposed here 
10000 loops, best of 3: 100 µs per loop