There are a small number of different Bs and different Ds, and I have to compute the quadratic form for a number of different combinations. I can probably do things the way you suggested, but I would lose the nice symmetry that I have in my current R code. (I am trying to switch to Julia for speed.)
On Monday, April 27, 2015 at 2:30:00 PM UTC-5, Andreas Noack wrote: > > If B and D are different then why is it not okay to calculate x = C\D and > then B'x afterwards? > > 2015-04-27 15:24 GMT-04:00 matt <[email protected] <javascript:>>: > >> I would like to compute multiple quadratic forms B'*C^(-1)*D, where B, C, >> and D are sparse, and C is always the same symmetric positive matrix. For >> fast computation, I would like to do this by computing and re-using the >> Cholesky factor: C=LL'. However, chol(C) does not work (because C is >> sparse), and cholfact(C)\D returns C^(-1)*D instead of C^(-1/2)*D. >> >> Minimal working example: >> C=4*speye(3); # sparse matrix >> B=sparse(ones(1,3)) # sparse unit vector >> chol(C) # returns "ERROR: Use cholfact()instead of chol() for sparse >> matrices..." >> L=cholfact(C); # works, but returns Cholesky decomposition instead of >> cholesky factor >> L\B' # would like this to be vector with elements .5 instead of .25 >> >> Thank you! >> > >
