Wicher2
commented
8 years ago Hi Haziq,
- Yes, this is a major bottleneck, especially for larger data sets, so these should be stored. This requires extra storage space, but that doesn't tend to be a problem in my experience. (For simple low-dimensional linear models this is not so eficient, but I wouldn't worry about that initially.)
- After 1. has been done, this becomes the main bottleneck I think. There is also a possible problem with numerical instability (whether using solve or choleski). In fact, the inverse does not need to be computed, since it is always multiplied by a vector anyway, which entails solving a linear equation, which is faster. There are two issues, computing the determinant of
Var.Y, and computingsolve(Var.Y) %*% afor a vectora.
For solving linear equations stably and fast (though I did initially make an error here as we found last week):
(\*solve (A+s\*I)x=b; stable for small s and singular A\*)
Clear[linsolve]
linsolve[a_, s_, b_] := Module[{eval, evec},
{eval, evec} = Eigensystem[a];
evec\[Transpose].((evec.b)/(eval + s))
];The following generates a linear solve function, so the eigensystem only needs to be computed once if eg both solve(A,b) and solve(A,c) need to be computed, best to implement this function but maybe try out the above first.
(\*create function u such that u[c] solves (A+s*I)x=c\*)
linsolve[a_, s_] := Module[{eval, evec},
{eval, evec} = Eigensystem[a];
evec\[Transpose].((evec.#)/(eval + s)) &
];Here & after an expression denotes a "pure function" with argument #.
Computing a determinant directly using in-built functions is potentially especially unstable (and I think both R and Mathematica use Intel libraries for this, so no fundamental differences). For the log determinant, the following is a standard relatively stable and fast method:
logdet[m_]:=2Total[Log[Diagonal[CholeskyDecomposition[(m+m\[Transpose])/2]]]]However, more efficient would be to write a procedure that combines the above linsolve and also compute the log determinant, in the above procedure a line like this would need to be added:
logdet=Total[Log[eval+s]]and give this as output as well. Still need to implement this myself.
Actually now did this and the following seems most convenient for the present purposes, so output the eigenvalues of H \Psi H (these can be used eg for computing determinants but also for Tr[Inverse[Var.Y]]) and the linearsolve function:
EvalAndLinsolve[a_, s_] := Module[{eval, evec},
{eval, evec} = Eigensystem[a];
{eval, evec\[Transpose].((evec.#)/(eval + s)) &}
];- 3 I always use
intercept=mean(Y), why might this not work? - 4 Yes, also here built in functions are slow and unstable for our setting. I think best to use the procedures under 2.
haziqj
After some timing tests, these are the bottlenecks in the programme:
H.mat %*% H.mat. Need some clever way to reduce these occurences. This is doable I think. The plan is to do all possible matrix multiplications and store them, and use them as needed. This way they are not evaluated again and again.Var.Y. The default R method is to usesolve(Var.Y). A faster, albeit less stable, way is to use the Cholesky decompositionchol2inv(chol(Var.Y)). This is 5 times faster than the regular method.Var.Y. If I can find another method which does not depend on this inverse, then I can speed up the programme. For example, sometimes it is as simple asintercept = mean(Y).dmvnorm()also takes a substantial amount of time, but not sure how much can be reduced by building own function. For reference,dmvnorm()takes 1/8th of the time to do asolve()Will look into this.@Wicher2 can you comment on the above points, specifically 2-4? Thanks.