# Most Efficient Way to Compute a Quadratic Matrix Form

**URL:** <https://discourse.julialang.org/t/most-efficient-way-to-compute-a-quadratic-matrix-form/66606>\
**Category:** Performance\
**Tags:** performance, linearalgebra\
**Created:** [August 18, 2021, 3:04pm UTC](https://discourse.julialang.org/t/most-efficient-way-to-compute-a-quadratic-matrix-form/66606 "2021-08-18T15:04:45Z")\
**Posts on this page:** 16\
**Page:** 1

<div class="post-metadata">

**Author:** ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)\
**Post date:** [August 18, 2021, 3:04pm UTC](https://discourse.julialang.org/t/most-efficient-way-to-compute-a-quadratic-matrix-form/66606/1 "2021-08-18T15:04:45Z")

</div>

Given a Symmetric Positive Semi definite Matrix P , what would be the most efficient (Fast) way to calculate the matrix quadratic form:

{x}^{T} P x

For the following cases:

1. The matrix P is dense.
2. The matrix P is sparse.
3. The matrix P is only symmetric (No definition).

---

<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:** [August 18, 2021, 3:34pm UTC](https://discourse.julialang.org/t/most-efficient-way-to-compute-a-quadratic-matrix-form/66606/2 "2021-08-18T15:34:30Z")

</div>

I think the answer will depend on how rank-deficient the matrix is (and its size), at least in the dense case. Is P reasonably close to full rank, or do you expect it to be severely degenerate?

~~If it is full rank, then I would claim that the fastest thing to do for (1) is `sum(abs2, cholesky(P).L'*x)`~~. And if the matrix is low enough rank, then I’d be tempted to compute a very high precision low rank approximation and then put that in. But again, there’s obviously overhead with something like that and so things like the size of P and how degenerate it is will presumably be important.

EDIT: sorry, I always compute x^T P^{-1} x, not the multiplication, so my solution for (1) was an autopilot thing and not going to be best for multiplications. And I even forgot to put the transpose in. Not my most precise comment, sorry.

---

<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:** [August 18, 2021, 5:40pm UTC](https://discourse.julialang.org/t/most-efficient-way-to-compute-a-quadratic-matrix-form/66606/3 "2021-08-18T17:40:09Z")

</div>

I would just use `dot(x, P, x)` (to compute x^\* P x; I’m not sure if you care whether it is conjugated in the complex case?), and worry about further micro-optimization only if benchmarks/profiling justifies the effort. This should be reasonably efficient in all cases even if it is not necessarily “optimal”.

---

<div class="post-metadata">

**Author:** ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)\
**Post date:** [August 18, 2021, 5:58pm UTC](https://discourse.julialang.org/t/most-efficient-way-to-compute-a-quadratic-matrix-form/66606/4 "2021-08-18T17:58:33Z")

</div>

You may assume it is something I found I need to optimize as much as possible.  
As far as I know, the `dot()` won’t take advantage of `mP` being a symmetric matrix.  
Is there an implementation which at least do that and gain some speed?

---

<div class="post-metadata">

**Author:** ![fipelle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fipelle/32/4772_2.png) [@fipelle](https://discourse.julialang.org/u/fipelle)\
**Post date:** [August 18, 2021, 6:19pm UTC](https://discourse.julialang.org/t/most-efficient-way-to-compute-a-quadratic-matrix-form/66606/5 "2021-08-18T18:19:22Z")

</div>

I have asked something similar recently ([Efficient approach to multiply three matrices (M1\*M2\*M3) and two vectors and a matrix (x\*M\*y)](https://discourse.julialang.org/t/efficient-approach-to-multiply-three-matrices-m1-m2-m3-and-two-vectors-and-a-matrix-x-m-y/66130)).

---

<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:** [August 18, 2021, 6:26pm UTC](https://discourse.julialang.org/t/most-efficient-way-to-compute-a-quadratic-matrix-form/66606/6 "2021-08-18T18:26:55Z")

</div>

Here’s an example where a matrix that is sufficiently large and rank-deficient benefits from something fancier:

```julia
using KernelMatrices, BenchmarkTools

const M = randn(50, 10_000);
const A = M'M;
const x = randn(10_000)

function fact_method(A,x)
  (U,V) = KernelMatrices.ACA(A, 1e-14)
  dot(U'*x, V'*x)
end

@btime dot($x, $A, $x) # 27.1 ms on my machine w/ 4 LAPACK threads
@btime fact_method($A, $x) # 27.2 ms with same setup.

```

I picked a rank that is about the break-even speed, but if you make it lower than 50 then the `fact_method` wins.

This factorization is the [adaptive cross approximation](https://link.springer.com/article/10.1007/s00607-002-1469-6), which I think is reasonably comparable to an LU factorization with a fancy pivoting strategy, but which chooses the pivots based on partial information to avoid passing over the whole matrix. I’m not super knowledgeable about the lowest level details about pivoting strategies, though, so that might be incorrect. In any case, I implemented it [here](https://bitbucket.org/cgeoga/kernelmatrices.jl/src/a0ce415f31e3f63aff76633056641ae93fca2cbf/src/factorizations.jl#lines-9) a long time ago. I don’t think it’s embarrassingly slow, but that package was one of my first big Julia projects and so while I haven’t seen anything better I bet somebody like @stevengj could sit down and write something much faster in the timespan of a coffee break.

I bring all this up to say that the code isn’t great and the factorization as currently implemented doesn’t specialize for symmetry, but even still for non-ridiculous matrix conditions it does beat the naive `dot(x, P, x)`. I would bet that a smarter implementation could be even more competitive. So I do think that the answer will really depend on some parameters of the problem like size and rank.

---

<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:** [August 18, 2021, 7:07pm UTC](https://discourse.julialang.org/t/most-efficient-way-to-compute-a-quadratic-matrix-form/66606/7 "2021-08-18T19:07:00Z")

</div>

> [@cgeoga](#):
>
> ```julia
> const M = randn(50, 10_000);
> const A = M'M;
> 
> ```

In this kind of problem, you’d be much better off not forming `A` explicitly in the first place…

---

<div class="post-metadata">

**Author:** ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)\
**Post date:** [August 18, 2021, 7:08pm UTC](https://discourse.julialang.org/t/most-efficient-way-to-compute-a-quadratic-matrix-form/66606/8 "2021-08-18T19:08:18Z")

</div>

I agree. If I had its factorization (Root, Cholesky, etc… or any other) it would have been easy.  
But assume I don’t have it.

---

<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:** [August 18, 2021, 7:18pm UTC](https://discourse.julialang.org/t/most-efficient-way-to-compute-a-quadratic-matrix-form/66606/9 "2021-08-18T19:18:39Z")

</div>

Obviously I’m not suggesting that you form the large rank-deficient square matrix and then re-factorize it. My point in that post is that the cost of the factorization might be worth it in some circumstances if the matrix is large enough rank deficient enough, and the matrix `A` is designed to be an example of that.

@RoyiAvital, can you comment on any more specifics for your problem, like the size and rank? I think without more specific information it’s hard to give particularly productive advice.

---

<div class="post-metadata">

**Author:** ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)\
**Post date:** [August 18, 2021, 7:38pm UTC](https://discourse.julialang.org/t/most-efficient-way-to-compute-a-quadratic-matrix-form/66606/10 "2021-08-18T19:38:47Z")

</div>

I think your direction is very well understood, if the rank is low enough one could apply factorization to speed things up (Cholesky for P being Positive Definite, LDL in case it is only symmetric would be classic choice while [`KernelMatrices.jl`](https://bitbucket.org/cgeoga/kernelmatrices.jl) may offer more advanced methods [See [`KernelMatrices.jl` - A Software Package for Working with Hierarchical Matrices](https://discourse.julialang.org/t/21696)]).

What about the case the rank is high to make those method not feasible.  
Is there any function which can take advantage of the structure of the problem?

---

<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:** [August 18, 2021, 8:02pm UTC](https://discourse.julialang.org/t/most-efficient-way-to-compute-a-quadratic-matrix-form/66606/11 "2021-08-18T20:02:36Z")

</div>

I appreciate the shout-out, although it’s really just the ACA that would be relevant for forming efficient low-rank representations like A = U V^T, and I have a feeling most people here could write a better ACA than what’s in there. So I would really only suggest anybody use that functionality as a cheap way to test for a proof of concept or something.

In any case, I don’t really have any ideas or knowledge about more thoughtful methods than the three-argument `dot`, at least at this level of generality and in a setting where factorizations aren’t helpful. But I’ll definitely be watching this thread, because I also find the topic interesting and post here in part because I’d love to see faster methods than the ones I currently use…

---

<div class="post-metadata">

**Author:** ![EricJohnson](https://avatars.discourse-cdn.com/v4/letter/e/ecae2f/32.png) [@EricJohnson](https://discourse.julialang.org/u/EricJohnson)\
**Post date:** [January 10, 2022, 2:33pm UTC](https://discourse.julialang.org/t/most-efficient-way-to-compute-a-quadratic-matrix-form/66606/12 "2022-01-10T14:33:19Z")

</div>

It seems there is a native implementation inside the standard library which is quite fast.  
At least according to the answer for [An Efficient Way (Fastest) to Calculate a Quadratic Form in Julia](https://stackoverflow.com/questions/70645845/an-efficient-way-fastest-to-calculate-a-quadratic-form-in-julia).

---

<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:** [January 10, 2022, 3:18pm UTC](https://discourse.julialang.org/t/most-efficient-way-to-compute-a-quadratic-matrix-form/66606/13 "2022-01-10T15:18:26Z")

</div>

For convenience I place it here (after some relabelling)

```julia
julia> using LinearAlgebra
julia> using BenchmarkTools

julia> x = rand(10000);
julia> A = rand(10000, 10000);

julia> P = (A + A') / 2;
julia> Pₛ = Symmetric(P);

julia> @btime dot($x, $P, $x);
  51.301 ms (0 allocations: 0 bytes)

julia> @btime dot($x, $Pₛ, $x);
  29.362 ms (0 allocations: 0 bytes)

```

---

<div class="post-metadata">

**Author:** ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)\
**Post date:** [January 11, 2022, 7:11am UTC](https://discourse.julialang.org/t/most-efficient-way-to-compute-a-quadratic-matrix-form/66606/14 "2022-01-11T07:11:53Z")

</div>

So it seems things are in good shape for the Dense case.  
Moreover, it seems that @Oscar_Smith found a way to [make things even faster](https://stackoverflow.com/questions/70645845/how-can-i-efficiently-calculate-a-quadratic-form-in-julia#comment124901999_70653944).

What about the sparse case? What does it use?

---

<div class="post-metadata">

**Author:** ![JM\_Beckers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jm_beckers/32/22482_2.png) [@JM\_Beckers](https://discourse.julialang.org/u/JM_Beckers)\
**Post date:** [January 11, 2022, 11:45am UTC](https://discourse.julialang.org/t/most-efficient-way-to-compute-a-quadratic-matrix-form/66606/15 "2022-01-11T11:45:53Z")

</div>

I think it also depends if you need to do the calculation only once, or if you need to repeat it for a lot of different vectors x but fixed P. In the latter case paying the price to calculate a reduced rank approximation could be worth doing.

---

<div class="post-metadata">

**Author:** ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)\
**Post date:** [October 27, 2024, 6:24pm UTC](https://discourse.julialang.org/t/most-efficient-way-to-compute-a-quadratic-matrix-form/66606/16 "2024-10-27T18:24:33Z")

</div>

I tried to replicate @Oscar_Smith 's result from the [StackOverflow comment](https://stackoverflow.com/questions/70645845/how-can-i-efficiently-calculate-a-quadratic-form-in-julia#comment124901999_70653944):

```julia
using BenchmarkTools;
using LinearAlgebra;
using LoopVectorization;

function QuadForm( vX :: Vector{T}, mA :: Matrix{T}, vY :: Vector{T} ) where {T <: AbstractFloat}

    (axes(vX)..., axes(vY)...) == axes(mA) || throw(DimensionMismatch());
    m, n = size(mA);
    s = zero(T);
    @tturbo for jj in 1:n
        yj = vY[jj];
        t = zero(T);
        for ii in 1:m
            t += mA[ii, jj] * vX[ii];
        end
        s += t * yj;
    end

    return s;

end

function QuadForm( vX :: Vector{T}, mA :: S, vY :: Vector{T} ) where {T <: AbstractFloat, S <: Symmetric{<: T, <: Matrix{<: T}}}
    # Slower than not using the Symmetric property

    (length(vX) == length(vY) == size(mA, 1)) || throw(DimensionMismatch())
    n = length(vY);
    s = zero(T);
    if mA.uplo == 'U'
        @inbounds for jj in 1:n
            @fastmath s += vX[jj] * mA[jj, jj] * vY[jj]; #<! Diagonal
            @inbounds for ii in 1:(jj - 1)
                @fastmath s += vX[ii] * mA[ii, jj] * vY[jj] + vX[jj] * mA[ii, jj] * vY[ii];
            end
        end
    else #<! if A.uplo == 'L'
        @inbounds for jj in 1:n
            @fastmath s += vX[jj] * mA[jj, jj] * vY[jj]; #<! Diagonal
            @inbounds for ii in (jj + 1):n
                @fastmath s += vX[ii] * mA[ii, jj] * vY[jj] + vX[jj] * mA[ii, jj] * vY[ii];
            end
        end
    end

    return s;

end

numRows = 10; numCols = 5; mA = randn(numRows, numCols); mB = Symmetric(randn(numRows, numRows)); mC = Matrix(mB); vX = randn(numRows); vY = randn(numCols);
numRows = 100; numCols = 50; mA = randn(numRows, numCols); mB = Symmetric(randn(numRows, numRows)); mC = Matrix(mB); vX = randn(numRows); vY = randn(numCols);
numRows = 1_000; numCols = 500; mA = randn(numRows, numCols); mB = Symmetric(randn(numRows, numRows)); mC = Matrix(mB); vX = randn(numRows); vY = randn(numCols);

@btime dot($vX, $mA, $vY)
@btime QuadForm($vX, $mA, $vY)
@btime dot($vX, $mB, $vX)
@btime QuadForm($vX, $mB, $vX) #<! Slow!
@btime dot($vX, $mC, $vX)
@btime QuadForm($vX, $mC, $vX) #<! Faster (No symmetry)

```

I’m on Ryzen 7940.  
On my system the Symmetric code make no gains.  
Up to 500 elements the `LoopVectorization.jl` code is 30-50% faster.

Is there anything else to optimize? Specifically for the symmetric case.  
Since the inner loop depends on the outer loop, `LP.jl` is not valid.