The frequency of people having these nearly symmetric (nearly Hermitian, etc.) matrices that they generated suggest to me that it would be nice to be able to apply the Symmetric wrapper to force symmetry (already works), but that we may want to consider the true symmetrized value to be the average of the two sides or something like that – after all, we have no real reason to choose the top or bottom half as correct.
On Wed, Nov 12, 2014 at 9:19 PM, Andreas Noack <[email protected] > wrote: > It is actually a bit tricky. Julia doesn't know that the matrix is > symmetric and she doesn't check for symmetry to save the time it takes to > perform the check (but she checks if the matrix is triangular). Because she > doesn't know that the matrix is symmetric, the inverse is calculated by a > general factorization (the LU) and here the exact symmetry is lost because > of small floating point errors. The check for positive definiteness is the > joint test that the matrix is symmetric (or Hermitian) and can be > factorized by the Cholesky factorization and therefore it returns false for > the inverse. > > One solution is to tell Julia that you know that the matrix is positive > definite before you calculate the inverse. Then she can exploit that > information and maintain the symmetry, e.g. > > julia> A = [i==j?1:0.3 for i=1:5, j=1:5] > 5x5 Array{Union(Float64,Int64),2}: > 1 0.3 0.3 0.3 0.3 > 0.3 1 0.3 0.3 0.3 > 0.3 0.3 1 0.3 0.3 > 0.3 0.3 0.3 1 0.3 > 0.3 0.3 0.3 0.3 1 > > julia> isposdef(inv(cholfact(A))) > true > > 2014-11-12 20:43 GMT-05:00 Sibnarayan Guria <[email protected]>: > > v=[i==j?1:0.3 for i=1:5, j=1:5] >> isposdef(v) ### true >> isposdef(inv(v)) ### false >> >> >> ---------- >> why? >> > >
