A.T * diag (b) * A + C 형태의 희소 행렬을 만드는 가장 빠른 방법은?

Question

내부 점 방식을 사용하여 대형 스파 스 비선형 시스템을 해결하는 코드를 최적화하려고합니다. 업데이트 단계에서 헤 시안 행렬 H, 그래디언트 g을 계산 한 다음 d을 H * d = -g으로 풀어 새로운 검색 방향을 얻습니다.

헤 시안 매트릭스 형태의 대칭 삼중 대각 구조를 갖는다 : * DIAG (B) * (A)에서

을 + C

I는 상기 특정 기능에 line_profiler을 실행 한 질문 :

Line # Hits  Time Per Hit % Time Line Contents 
================================================== 
    386        def _direction(n, res, M, Hsig, scale_var, grad_lnprior, z, fac): 
    387        
    388         # gradient 
    389 44 1241715 28220.8 3.7  g = 2 * scale_var * res - grad_lnprior + z * np.dot(M.T, 1./n) 
    390        
    391         # hessian 
    392 44 3103117 70525.4 9.3  N = sparse.diags(1./n ** 2, 0, format=FMT, dtype=DTYPE) 
    393 44 18814307 427597.9 56.2  H = - Hsig - z * np.dot(M.T, np.dot(N, M)) # slow! 
    394         
    395         # update direction 
    396 44 10329556 234762.6 30.8  d, fac = my_solver(H, -g, fac) 
    397         
    398 44  111  2.5 0.0  return d, fac 

출력을 보면 H은 훨씬 비용이 많이 드는 단계입니다. 실제로는 새로운 방향을 위해 실제로 해결하는 것보다 훨씬 오래 걸립니다.

Hsig 및 M 모두 CSC 스파 스 매트릭스이며, n는 밀도 벡터이며 z는 스칼라이다. 사용하고있는 해답은 H이 CSC 또는 CSR 스파 스 매트릭스가되어야합니다. 여기

는 내 진짜 매트릭스와 같은 형식, 크기 및 희소과 일부 장난감 데이터를 생성하는 기능입니다 :

def original(Hsig, M, n, z): 
    N = sparse.diags(1./n ** 2, 0, format='csc') 
    H = - Hsig - z * np.dot(M.T, np.dot(N, M)) # slow! 
    return H 

타이밍 : 여기

import numpy as np 
from scipy import sparse 

def make_toy_data(nt=200000, nc=10): 

    d0 = np.random.randn(nc * (nt - 1)) 
    d1 = np.random.randn(nc * (nt - 1)) 
    M = sparse.diags((d0, d1), (0, nc), shape=(nc * (nt - 1), nc * nt), 
        format='csc', dtype=np.float64) 

    d0 = np.random.randn(nc * nt) 
    Hsig = sparse.diags(d0, 0, shape=(nc * nt, nc * nt), format='csc', 
         dtype=np.float64) 

    n = np.random.randn(nc * (nt - 1)) 
    z = np.random.randn() 

    return Hsig, M, n, z 

그리고 H을 구성하는 내 원래의 접근 방법이있다 :

%timeit original(Hsig, M, n, z) 
# 1 loops, best of 3: 483 ms per loop 

빠른 방법이 있습니까 이 행렬?

Answer 1

입니다.d0 및 d1이 M의 주요 상부 대각선, 그리고 d이 D의 주요 대각선 경우, 다음 코드는 M.T * D * M 직접 작성하여 매트릭스 M는 CSR 형식이었다

def make_tridi_bis(d0, d1, d, nc=10): 
    d00 = d0*d0*d 
    d11 = d1*d1*d 
    d01 = d0*d1*d 
    len_ = d0.size 
    data = np.empty((3*len_ + nc,)) 
    indices = np.empty((3*len_ + nc,), dtype=np.int) 
    # Fill main diagonal 
    data[:2*nc:2] = d00[:nc] 
    indices[:2*nc:2] = np.arange(nc) 
    data[2*nc+1:-2*nc:3] = d00[nc:] + d11[:-nc] 
    indices[2*nc+1:-2*nc:3] = np.arange(nc, len_) 
    data[-2*nc+1::2] = d11[-nc:] 
    indices[-2*nc+1::2] = np.arange(len_, len_ + nc) 
    # Fill top diagonal 
    data[1:2*nc:2] = d01[:nc] 
    indices[1:2*nc:2] = np.arange(nc, 2*nc) 
    data[2*nc+2:-2*nc:3] = d01[nc:] 
    indices[2*nc+2:-2*nc:3] = np.arange(2*nc, len_+nc) 
    # Fill bottom diagonal 
    data[2*nc:-2*nc:3] = d01[:-nc] 
    indices[2*nc:-2*nc:3] = np.arange(len_ - nc) 
    data[-2*nc::2] = d01[-nc:] 
    indices[-2*nc::2] = np.arange(len_ - nc ,len_) 

    indptr = np.empty((len_ + nc + 1,), dtype=np.int) 
    indptr[0] = 0 
    indptr[1:nc+1] = 2 
    indptr[nc+1:len_+1] = 3 
    indptr[-nc:] = 2 
    np.cumsum(indptr, out=indptr) 

    return sparse.csr_matrix((data, indices, indptr), shape=(len_+nc, len_+nc)) 

경우 d0을 추출 할 수 있습니다 및 d1d0 = M.data[::2] 및 d1 = M.data[1::2], 나는 내가 무엇을 얻을 당신 장난감 데이터뿐만 아니라 그 배열을 반환하는 루틴을 만들고, 여기 있어요 수정 :

In [90]: np.allclose((M.T * sparse.diags(d, 0) * M).A, make_tridi_bis(d0, d1, d).A) 
Out[90]: True 

In [92]: %timeit make_tridi_bis(d0, d1, d) 
10 loops, best of 3: 124 ms per loop 

In [93]: %timeit M.T * sparse.diags(d, 0) * M 
1 loops, best of 3: 501 ms per loop 

---

위 코드의 전체 목적은 0이 아닌 항목의 구조를 이용하는 것입니다. 당신이 함께 증식하는 행렬의 다이어그램을 그릴 경우, (d_0) 메인 상단과 하단 결과 삼중 대각 행렬의 (d_1) 대각선 단순히 것을 자신을 설득하기가 상대적으로 쉽다 :

d_0 = np.zeros((len_ + nc,)) 
d_0[:len_] = d00 
d_0[-len_:] += d11 

d_1 = d01 

그 함수의 나머지 코드는 위의 데이터를 가진 sparse.diags을 호출하는 것이 몇 배 더 느리기 때문에 단순히 삼중 대 매트릭스를 직접 구축합니다.

Answer 2

테스트 케이스를 실행 해 보았는데 np.dot(N, M)에 문제가있었습니다. 나는 그걸 파고 들지 않았지만, 내 numpy/sparse 콤보 (둘다 꽤 새로운)는 희소 배열에 np.dot을 사용하는 데 문제가 있다고 생각합니다.

그러나 H = -Hsig - z*M.T.dot(N.dot(M))은 정상적으로 실행됩니다. 여기에는 sparse dot이 사용됩니다.

프로파일을 실행하지 않았지만 여기에 여러 부분에 대한 Ipython 타이밍이 있습니다. 그 이중 점을 만드는 것보다 데이터를 생성하는 데 시간이 오래 걸립니다.

In [37]: timeit Hsig,M,n,z=make_toy_data() 
1 loops, best of 3: 2 s per loop 

In [38]: timeit N = sparse.diags(1./n ** 2, 0, format='csc') 
1 loops, best of 3: 377 ms per loop 

In [39]: timeit H = -Hsig - z*M.T.dot(N.dot(M)) 
1 loops, best of 3: 1.55 s per loop 

H 내가 세 대각선 배열에서 제품 M.T * D * M을 계산에 4 배의 속도 향상에 가까운 얻을

<2000000x2000000 sparse matrix of type '<type 'numpy.float64'>' 
    with 5999980 stored elements in Compressed Sparse Column format>