# Representing a symmetric matrix factorization
**URL:** https://discourse.julialang.org/t/representing-a-symmetric-matrix-factorization/122881
**Category:** Numerics
**Tags:** question, linearalgebra
**Created:** [November 21, 2024, 8:56am UTC](https://discourse.julialang.org/t/representing-a-symmetric-matrix-factorization/122881 "2024-11-21T08:56:58Z")
**Posts on this page:** 4
**Page:** 1
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [November 21, 2024, 8:56am UTC](https://discourse.julialang.org/t/representing-a-symmetric-matrix-factorization/122881/1 "2024-11-21T08:56:58Z")
This question is about implementing a matrix factorization with some operations in Julia, ideally relying on existing optimized implementations, in a generic way do that I don’t have to commit to a particular implementation (dense/sparse, etc) and have it work with AD.
### What I need
Suppose J is an n \times m matrix (n \ge m possible). We form B = J'J, and want to support the following operations:
1. products x'Bx for any conformable x vectors,
2. testing if B is positive definite (it is always positive semidefinite by construction),
3. B^{-1} y for conformable y vectors.
4. Forming the actual matrix B (for input into `hvcat`).
### Implementation questions
Ideally these should be fast, work for dense and sparse B (J is a Jacobian obtained from AD, may be sparse). The code should be generic and work with AD (`ForwardDiff.Dual`, Enzyme.jl, etc).
I thought of making a J=QR decomposition and storing R. This is not unlike a Cholesky decomposition, except that `LinearAlgebra.Cholesky` requires positive definiteness if I understand correctly, and `dot(x, B, x)` does not specialize for `Cholesky` (but the generic fallback may not be so bad). And I did not check AD’ability of the sparse representation.
I am wondering if there is a better way to do it though — can one arrive at an LDL' representation from J without forming B? Or is there an better factorization?
---
### Author:  [@tim.holy](https://discourse.julialang.org/u/tim.holy)
#### Post date: [November 21, 2024, 9:10am UTC](https://discourse.julialang.org/t/representing-a-symmetric-matrix-factorization/122881/2 "2024-11-21T09:10:27Z")
Sounds like you might want `ldlt`? With a wrapper that’s a lot like [`PDMat`](https://github.com/JuliaStats/PDMats.jl/blob/592634b4dfadcf9ed74a6123ceff23c776bd93a1/src/pdmat.jl#L1-L9) except with mathematics appropriate for Bunch-Kauffman?
> (it is always positive semidefinite by construction)
Though if this is true, why not just use `PDMat`?
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [November 21, 2024, 9:30am UTC](https://discourse.julialang.org/t/representing-a-symmetric-matrix-factorization/122881/3 "2024-11-21T09:30:32Z")
> [@tim.holy](#):
>
> Sounds like you might want `ldlt`?
Possibly, but recall that the input is J, and I am not sure if I want to form B = J'J. Is that the only way to get B = LDL'?
> [@tim.holy](#):
>
> Though if this is true, why not just use `PDMat`?
Because it forms both the decomposition and the full matrix unconditionally. I rarely need the full matrix.
---
### Author:  [@tim.holy](https://discourse.julialang.org/u/tim.holy)
#### Post date: [November 21, 2024, 9:56am UTC](https://discourse.julialang.org/t/representing-a-symmetric-matrix-factorization/122881/4 "2024-11-21T09:56:35Z")
Ah, I didn’t understand that `J` was supplied. I agree with computing the QR decomposition and storing `R`. Then
```julia
struct SymProduct{T,RT<:AbstractMatrix{T}} <: Factorization{T}
R::RT
end
Base.size(B::SymProduct) = (m = size(B.R, 1); return (m, m))
Base.:*(B::SymProduct, x) = (y = B.R*x; return y'*y)
function LinearAlgebra.isposdef(B::SymProduct)
(; R = B)
m, n = size(R)
m > n && return false # fast path
# check for zeros in R
end
...
```
If you want `SymProduct{T,RT<:AbstractMatrix{T}} <: AbstractMatrix{T}` instead, then you need to store `B` because `B[i, j]` needs to be fast for it to truly be considered an `AbstractMatrix` (that’s basically the defining difference between `Factorization` and `AbstractMatrix`).
It would be useful to have this as a standalone package, as this does come up in multiple applications. JuliaLinearAlgebra could host, ping me here or on slack if you need permissions. Or perhaps it could be added to PDMats.jl.