Package for dealing with symmetric matrices efficiently

The Julia Symmetric Matrix type stores the full matrix. I’m looking for a package that offers a way to store and to index into a symmetric matrix directly but only stores one triangle, but that is performant, operates as a matrix, and also allows me to change off-diagonal elements at will, e.g. setting X[2,3] = 4 is an acceptable operation.

What is my best bet here?

Background: storing Hessians for large-dimensional optimization problems.

Thanks!!

3 Likes

Google serves up StandardPacked.jl and RectangularFullPacked.jl. Note, as mentioned in the READMEs, that the memory savings come at the cost of some extra overhead and a slight loss of floating-point accuracy in many calculations.

3 Likes

Thanks a bunch @danielwe . I’ll have a look at these too. There’s a whole bunch of them, many dormant, now incompatible, or restrictive.

@danielwe Do you have a recommendation?

Nope, never used any of these.

Thanks. Guess I should have phrased my original question differently. I’m looking for which package to use out of the umpteen different ones that claim to do this.

I think the conventional wisdom here is that the factor-of-two savings of a packed format aren’t usually worth the performance hit that this incurs. Is a factor of two in memory really the difference between solvable and unsolvable for you?

Have you considered ways to avoid explicit storage of the whole n \times n Hessian? Note that you can compute Hessian–vector products much more efficiently without storing the matrix — in fact, you can save a factor of roughly n by forward-over-reverse-mode computation of Hessian–vector products (see section 8.4.1 of our matrix-calculus course notes)!

Given a fast Hessian–vector product, you can construct a low-rank approximation of the Hessian using e.g. randomized sketching methods in LowRankApprox.jl. (And L-BFGS optimization algorithms estimate low-rank Hessian approximations by a different method.)

6 Likes

In another direction, is it possible that your Hessian matrix is sparse?

1 Like

Thanks @stevengj . That’s helpful. No, it isn’t: waste not want not and all. :wink: In fact, storing is what I’ll do since I suspect that they are not large by your standards.

I’m aware of the LBFGS algorithms with low-rank Hessians. I’m using Newton with Trust Regions which appears to be better behaved for the type of problem I’m working on.

No, but thanks.