[grajlord]

Próbuję zaimplementować w R kod z MATLABa wykonujący transformację na wektorze. Metoda opisana jest tutaj - usuwa z wektora (szeregu czasowego; zmienna z w kodzie) trend zgodnie z zadanym parametrem lambda.

W sumie to mam problem z matematyką, a nie implamentacją. Jeżeli wektor z ma T długości (a z reguły będzie długi -
T > 100000), to macierze w obliczeniach będą wielkości zbliżonej do TxT, czyli ogromne. Macierz I jest na szczęście jednostkowa, a macierze D2 i D2' niezerowe tylko na odpowiednio 3 przekątnych, a do tego mam metody, które zjedzą niewiele pamięci. Wszystko się natomiast sypie, gdy trzeba znaleźć macierz odwrotną: (I + lambda^2 * D2' x D2)^-1, i to zarówno jeśli chodzi o pamięć, jak i prędkość obliczeń.
Czy jest jakaś właściwość macierzy, która pozwoli to sprawnie policzyć, albo ewentualnie jakiś sposób na efektywne przybliżenie tego wyniku? Tak na ludzki rozum, to skoro w szeregu czasowym największa korelacja istnieje między najbliższymi pomiarami, jednoczesne analizowanie 1-ego i T-ego punktu niewiele zmieni w wyniku.
Mógłbym podzielić z na mniejsze zachodzące na siebie fragmenty, policzyć to w pętli, a potem uśrednić 'ogony' z każdej pętli, ale może jest ładniejsze rozwiązanie?
Wiem, że metoda jest zaimplementowana w komercyjnym oprogramowaniu i działa szybko, dla T = 1e5 policzy w kilka sekund na domowym sprzęcie.

W matlabie idzie to tak:

T = length(z);
lambda = 10;
I = speye(T);
D2 = spdiags(ones(T-2,1)*[1 -2 1],[0:2],T-2,T);
z_stat = (I-inv(I+lambda^2*D2'*D2))*z;

W R byłoby tak:

require(Matrix)
z = rnorm(1e5, 1)
t = length(z)
lambda = 10
I = .sparseDiagonal(t)
mat <- matrix(1, t-2,1) %*% t(c(1, -2, 1))
D2 = bandSparse(t-2, t, 0:2, mat)
z_stat = (I - solve(I + lambda^2 * t(D2) %*% D2)) %*% z

pomocy @scyth i inni
#matematyka #r #matlab

[grajlord]

@grajlord: Dobra, już wiem. Jak gdzieś trzeba obliczać macierz odwrotną, to znaczy, żeby jej nie liczyć, tylko układ równać rozwiązać.

Czyli w R zamiast ostatniej linii

z_stat = (I - solve(I + lambda^2 * t(D2) %*% D2)) %*% z

wrzucam:

sparse = I + lambda^2 * t(D2) %*% D2
z_stat = z - solve(sparse, z)

W Matlabie nie wiem.

[scyth]

@grajlord: No ja nie wiem... :(

[grajlord]

@scyth: Na razie działa :P Trochę sobie potestowałem i wyszło na to, że przy zgodnych parametrach ta metoda daje wyniki bardzo zbliżone do loess. Z tą różnicą, że dla dużych setów liczy się zdecydowanie szybciej.