How large is your system? The calculation U * diagm(D) * U' * b could easily use all of your RAM if the problem is large. Even if it is not really large you should evaluate it from the right and use Diagonal, i.e. U * (Diagonal(D) * (U' * b)) because the you are not construction a large matrix, but only vectors.
If the problem is not that large, you might just want to write pinv(full(A))*b because that is effectively the same as what you are doing now just with a sparse solver. I don't think it is possible to specify a threshold in eigs. 2015-04-27 15:12 GMT-04:00 Mladen Kolar <[email protected]>: > Dear Andreas, > > Thank you for your response and I apologize for not writing earlier as I > got swamped with teaching. > > Your suggestion of using cholfact together with "shift" argument is a > reasonable one. I need to test it and see what results I'll get for a > specific problem at hand. Also, I should probably make a leap and install > v0.4. > > I would like to ask a related question on a slightly different approach. > While in general the solution is not unique when A is low rank, it seems > that in my problem vector b is in the range of A (this is an empirical > observation). The following approach seems to result in an estimator with > good statistical properties. > > # A -- sparse matrix > # b -- vector > > D, U, _ = eigs(A; nev = size(A, 1)-1 ) > indC = findfirst((x) -> x < 5e-5, D) > D[indC:end] = 0 > for j=1:indC > D[j] = 1. / D[j] > end > hat_x = U * diagm(D) * U' * b > > In the above snippet, first eigen-decomposition of A is performed. Next, a > low rank approximation of A is formed by keeping all eigenvalues larger > that 5e-5 (this is a tuning parameter). Finally, x is obtained by computing > the pseudo-inverse of the low-rank approximation and multiplying with b. > > The function eigs accepts the parameter "nev", however, in general I do > not know how many eigenvalues are going to be above a certain cutoff (in > the snippet above, 5e-5). Is it possible to pass a parameter to obtain all > eigenvalues and eigenvectors above certain threshold? Or should I simply > identify one eigenvalue at a time and deflate A by removing the > corresponding eigenvector? > > Thanks, Mladen > > > On Tuesday, April 21, 2015 at 9:09:29 AM UTC-5, Andreas Noack wrote: > >> I'm not sure what the best solution is here because I don't fully >> understand your objective. If A has low rank then the solution is not >> unique and if it has almost low rank, the solution is very ill conditioned. >> A solution could our new "shift" argument to our complex Cholesky >> factorization. This will introduce some regularization to the problem. This >> requires latest master of Julia. Then you can do >> julia> A = sprandn(20,20,0.1);A = A'A; >> >> julia> b = A*ones(size(A, 1)); >> >> julia> A\b >> ERROR: ArgumentError: matrix has one or more zero pivots >> in ldltfact at sparse/cholmod.jl:1042 >> in ldltfact at sparse/cholmod.jl:1048 >> in factorize at sparse/linalg.jl:688 >> in \ at linalg/generic.jl:244 >> >> julia> cholfact(A, shift = 1e-16)\b >> 20-element Array{Float64,1}: >> 1.0 >> 1.0 >> 1.0 >> 1.0 >> 1.0 >> 1.0 >> 1.0 >> 1.0 >> 1.0 >> 1.0 >> 1.0 >> 1.0 >> 1.38612 >> 0.0 >> 1.0 >> 1.0 >> 1.0 >> 1.0 >> 1.0 >> 1.0 >> >> 2015-04-20 20:32 GMT-04:00 Mladen Kolar <[email protected]>: >> >>> Hi, >>> >>> I am interested in solving a linear system Ax = b where A is a sparse, >>> symmetric positive semi-definite matrix. >>> A is not exactly low rank, however, for the estimation problem at hand, >>> one can approximate A with a matrix that has rank around sqrt(p) where A is >>> p x p matrix and still obtain an estimator of x that will have good >>> statistical properties. >>> >>> Is there an incomplete Cholesky factorization for sparse symmetric >>> matrices in julia? Or an svd implementation for sparse matrices? >>> What would be a recommended way of solving the above problem? >>> >>> Kind regards, >>> Mladen >>> >> >>
