# Fast computation of Grammian of matrix

**URL:** <https://discourse.julialang.org/t/fast-computation-of-grammian-of-matrix/56244>\
**Category:** Performance\
**Tags:** question\
**Created:** [March 1, 2021, 11:52am UTC](https://discourse.julialang.org/t/fast-computation-of-grammian-of-matrix/56244 "2021-03-01T11:52:55Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![BambOoxX](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bambooxx/32/22179_2.png) [@BambOoxX](https://discourse.julialang.org/u/BambOoxX)\
**Post date:** [March 1, 2021, 11:52am UTC](https://discourse.julialang.org/t/fast-computation-of-grammian-of-matrix/56244/1 "2021-03-01T11:52:56Z")

</div>

For some model indentification problems, I have to compute the grammian of data matrices i.e. compute `A = X'*X` where X is a complex rectangular matrix.

Is there a fast implementation of this specific operation in julia knowing that the output matrix is hermitian ?

A first basic assessment would assume this could split in half the computation times as well as the memory requirements. Also this could reduce numerical error propagation leading to a non heremitian numeric output.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [March 1, 2021, 12:00pm UTC](https://discourse.julialang.org/t/fast-computation-of-grammian-of-matrix/56244/2 "2021-03-01T12:00:45Z")

</div>

Try to search for “kernel matrices julia” and you may find something useful. For example, I found [KernelMatrices.jl](https://bitbucket.org/cgeoga/kernelmatrices.jl/src/master/) by @cgeoga

You may also find something if you search for “quadratic form julia”, I don’t know.

---

<div class="post-metadata">

**Author:** ![FMeirinhos](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fmeirinhos/32/19527_2.png) [@FMeirinhos](https://discourse.julialang.org/u/FMeirinhos)\
**Post date:** [March 1, 2021, 12:30pm UTC](https://discourse.julialang.org/t/fast-computation-of-grammian-of-matrix/56244/3 "2021-03-01T12:30:32Z")

</div>

Some quick googling finds the Julia benevolent overlords had already gotten you covered almost 10 years ago  
[https://github.com/JuliaLang/julia/issues/659](https://github.com/JuliaLang/julia/issues/659)

```julia
using BenchmarkTools
a = rand(100,100); b = rand(100,100);

@benchmark $a' * $a
# median time: 110.190 μs (0.00% GC)

@benchmark $b' * $a
# median time: 836.717 μs (0.00% GC)
```

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [March 1, 2021, 12:34pm UTC](https://discourse.julialang.org/t/fast-computation-of-grammian-of-matrix/56244/4 "2021-03-01T12:34:29Z")

</div>

Just a remark: it is often advantageous from a numerical viewpoint to use algorithms that do not require forming the gramian and instead work with its square root (Cholesky) factor.

---

<div class="post-metadata">

**Author:** ![BambOoxX](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bambooxx/32/22179_2.png) [@BambOoxX](https://discourse.julialang.org/u/BambOoxX)\
**Post date:** [March 1, 2021, 12:48pm UTC](https://discourse.julialang.org/t/fast-computation-of-grammian-of-matrix/56244/5 "2021-03-01T12:48:47Z")

</div>

@zdenek_hurak could you expand on your comment ? What type of algorithms do you refer to ?

---

<div class="post-metadata">

**Author:** ![BambOoxX](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bambooxx/32/22179_2.png) [@BambOoxX](https://discourse.julialang.org/u/BambOoxX)\
**Post date:** [March 1, 2021, 12:49pm UTC](https://discourse.julialang.org/t/fast-computation-of-grammian-of-matrix/56244/6 "2021-03-01T12:49:31Z")

</div>

@FMeirinhos thanks, sometimes it’s just a matter of keywords !

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [March 1, 2021, 1:02pm UTC](https://discourse.julialang.org/t/fast-computation-of-grammian-of-matrix/56244/7 "2021-03-01T13:02:46Z")

</div>

I meant something in the line of square-root balancing used for model order reduction ([Redirecting](https://doi.org/10.1016/S1474-6670(17)49168-3) or [Balancing free square-root algorithm for computing singular perturbation approximations | IEEE Conference Publication | IEEE Xplore](https://doi.org/10.1109/CDC.1991.261486) and certainly described elsewhere). Perhaps it may be relevant for you problem, maybe not. Just in case, the author has been actively developing [MatrixEquations.jl](https://github.com/andreasvarga/MatrixEquations.jl) package. Have a look at it, perhaps you may find something relevant there.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 1, 2021, 1:45pm UTC](https://discourse.julialang.org/t/fast-computation-of-grammian-of-matrix/56244/8 "2021-03-01T13:45:01Z")

</div>

> [@FMeirinhos](#):
>
> ```julia
> @benchmark $a' * $a
> # median time: 110.190 μs (0.00% GC)
> 
> @benchmark $b' * $a
> # median time: 836.717 μs (0.00% GC)
> 
> ```

Typically you should look at the minimum time rather than the median, since timing noise is always positive. (I usually use `@btime`, which only reports the minimum.) When you do that, the conclusion is the opposite on my computer:

```julia
julia> @btime $a' * $a;
  49.545 μs (2 allocations: 78.20 KiB)

julia> @btime $b' * $a;
  30.643 μs (2 allocations: 78.20 KiB)

```

Similarly for a larger matrix:

```julia
julia> @btime $A' * $A;
  342.402 μs (2 allocations: 78.20 KiB)

julia> @btime $B' * $A;
  246.269 μs (2 allocations: 78.20 KiB)

```

It looks like the culprit is threading — it seems that dgemm is better-parallelized than dysrk in OpenBLAS, probably because the computation in dgemm is more regular. If I change to single-threaded BLAS then dsyrk is faster in both cases:

```julia
julia> LinearAlgebra.BLAS.set_num_threads(1)

julia> @btime $A' * $A;
  487.042 μs (2 allocations: 78.20 KiB)

julia> @btime $B' * $A;
  602.739 μs (2 allocations: 78.20 KiB)

julia> @btime $a' * $a;
  58.816 μs (2 allocations: 78.20 KiB)

julia> @btime $b' * $a;
  61.733 μs (2 allocations: 78.20 KiB)

```

---

<div class="post-metadata">

**Author:** ![BambOoxX](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bambooxx/32/22179_2.png) [@BambOoxX](https://discourse.julialang.org/u/BambOoxX)\
**Post date:** [March 1, 2021, 2:02pm UTC](https://discourse.julialang.org/t/fast-computation-of-grammian-of-matrix/56244/9 "2021-03-01T14:02:38Z")

</div>

Thanks, I am not sure I fully understand the whole process yet, but it may be useful later.

---

<div class="post-metadata">

**Author:** ![cgeoga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cgeoga/32/216186_2.png) [@cgeoga](https://discourse.julialang.org/u/cgeoga)\
**Post date:** [March 1, 2021, 4:56pm UTC](https://discourse.julialang.org/t/fast-computation-of-grammian-of-matrix/56244/10 "2021-03-01T16:56:02Z")

</div>

To step back a little bit, do you actually need to compute that matrix? If you just need a few entries or to apply it to vectors or something, something like this pseudo-julia might be helpful:

```julia
struct GramMatrix{T} 
  X::Matrix{T} # Or perhaps a QR factorization of X?
end

Base.getindex(G::GramMatrix{T}, j::Int64, k::Int64) where{T} = # small inner product of j/k slices
Base.*(G::GramMatrix{T}, x::Vector{T}) where{T} = # manual two-part matrix multiplication
Base.\(G::GramMatrix{T}, x::Vector{T}) where{T} = ... #, perhaps using IterativeSolvers.jl with your fast matvec

# ....

```

Of course if you need the whole matrix then this won’t be enough, but particularly for very rectangular X you could potentially stand to save a lot with just a few lines of code like this if you only need something simple.

---

<div class="post-metadata">

**Author:** ![BambOoxX](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bambooxx/32/22179_2.png) [@BambOoxX](https://discourse.julialang.org/u/BambOoxX)\
**Post date:** [April 2, 2021, 3:31pm UTC](https://discourse.julialang.org/t/fast-computation-of-grammian-of-matrix/56244/11 "2021-04-02T15:31:15Z")

</div>

Sorry for the (very) late reply !

Actually, I am computing some reduced normal equations from covariance and cross covariance matrices.  
So I have some `CovX = X'*X`,`CrcovXY = X'*Y` and `CovY = Y'*Y` matrices, which are then assembled as  
`Z = CovY - CrcovXY'*CovX*CrcovY`.

Of course, this gets me in much numerical trouble and there seems to be a lot of error propagation.  
For instance, though `Z` should theoretically be Hermitian, it is numerically not.

---

<div class="post-metadata">

**Author:** ![BambOoxX](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bambooxx/32/22179_2.png) [@BambOoxX](https://discourse.julialang.org/u/BambOoxX)\
**Post date:** [April 6, 2021, 1:16pm UTC](https://discourse.julialang.org/t/fast-computation-of-grammian-of-matrix/56244/12 "2021-04-06T13:16:09Z")

</div>

I read the articles you mention a bit more thoroughly. It seems that the available algorithms to perform the computation of the Cholesky decomposition of the Gramian matrix require to have a realization of a system (whether in continuous or discrete time). Because of that, I do not see how to take advantage of this (yet promising otherwise) approach.