# Faster quadratic expression for symmetric matrices

**URL:** <https://discourse.julialang.org/t/faster-quadratic-expression-for-symmetric-matrices/91843>\
**Category:** Performance\
**Tags:** question, linearalgebra\
**Created:** [December 19, 2022, 11:05am UTC](https://discourse.julialang.org/t/faster-quadratic-expression-for-symmetric-matrices/91843 "2022-12-19T11:05:32Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![corbat](https://avatars.discourse-cdn.com/v4/letter/c/e0b2c6/32.png) [@corbat](https://discourse.julialang.org/u/corbat)\
**Post date:** [December 19, 2022, 11:05am UTC](https://discourse.julialang.org/t/faster-quadratic-expression-for-symmetric-matrices/91843/1 "2022-12-19T11:05:32Z")

</div>

Hi, I have a use case where I need to evaluate many quadratic expressions of the form x^T M x, with the matrix M always beeing symmetric. The typical dimensionality of M is between 3 and 3000.

Using the following two approaches for computing this expresssion,

```julia
product_1(x, M) = dot(x, M, x)

function product_2(x, M) 
	Mx = M * x
	return dot(x, Mx)
end 

```

and when comparing `Matrix{Float64}` and `Symmetric{Float64, Matrix{Float64}}` matrices, I found the following performances:

 ![btime_comparison](https://global.discourse-cdn.com/julialang/original/3X/5/5/550df5e552ae5417c34161676c451032c3e31cb1.png)

I’m a bit surprised that the `Symmetric` case is not always faster. Also, I’m not sure how “bad” it is to make the tradeoff between computation time and memory allocations when using `product_2` for large matrices (for D=1000: 1 allocation: ~ 8 KiB, compared to 0 allocations for `product_1`)

I also tried to implement this product myself, since due to the symmetry of M, only either the lower or the upper triangle of M needs to be considered. However, I have not been able to make this product as fast as the other two functions, even though I only iterate over half of the matrix.

```julia
function product_3(x, M)
    n = length(x)
    result = zero(eltype(x))

    @inbounds for i = 1:n
        result += x[i]^2 * M[i, i]

        @simd for j = 1 : i - 1
            result += 2* x[i] * x[j] * M[j, i]
        end 
    end

    return result
end

```

(`product_3` is much slower and has a lot of allocations.)

Do you have any suggestions on how to get the best performance for this type of product with symmetric matrices? Is there a way around using different functions and datatypes for different sizes of M to get the best performance?

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [December 19, 2022, 1:00pm UTC](https://discourse.julialang.org/t/faster-quadratic-expression-for-symmetric-matrices/91843/2 "2022-12-19T13:00:08Z")

</div>

it might be worth trying mkl

---

<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:** [December 19, 2022, 1:57pm UTC](https://discourse.julialang.org/t/faster-quadratic-expression-for-symmetric-matrices/91843/3 "2022-12-19T13:57:03Z")

</div>

> [@corbat](#):
>
> I also tried to implement this product myself, since due to the symmetry of M, only either the lower or the upper triangle of M needs to be considered.

The LinearAlgebra library already implements [exactly such a method](https://github.com/JuliaLang/julia/blob/427432e5c6ea90aa2f4616a380b4f4322ff30bbe/stdlib/LinearAlgebra/src/symmetric.jl#L579-L602).

There was some discussion of the performance issues when this was implemented, where it was noted that this was not as fast (for large matrices) as the naive approach using the BLAS matrix–vector multiplication. See [julia#32739 (add generalized dot product)](https://github.com/JuliaLang/julia/pull/32739).

---

<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:** [December 19, 2022, 2:21pm UTC](https://discourse.julialang.org/t/faster-quadratic-expression-for-symmetric-matrices/91843/4 "2022-12-19T14:21:50Z")

</div>

To do better here, you basically have to re-implement the BLAS `dsymv` (symmetrix matrix–vector product) function in `dot(x, A, b)`.

Optimizing this kind of thing is non-trivial, but to give you an idea, here is an optimized (but not fully debugged) _skew_-symmetric matrix–vector product in Julia developed by Simon Mataigne with help from @Elrod: [SkewLinearAlgebra.jl/skewsymmetric.jl at 54067de46fbe5dc3f1dd48546a64dc7b406b48b5 · JuliaLinearAlgebra/SkewLinearAlgebra.jl · GitHub](https://github.com/JuliaLinearAlgebra/SkewLinearAlgebra.jl/blob/54067de46fbe5dc3f1dd48546a64dc7b406b48b5/src/skewsymmetric.jl)

---

<div class="post-metadata">

**Author:** ![corbat](https://avatars.discourse-cdn.com/v4/letter/c/e0b2c6/32.png) [@corbat](https://discourse.julialang.org/u/corbat)\
**Post date:** [December 19, 2022, 3:49pm UTC](https://discourse.julialang.org/t/faster-quadratic-expression-for-symmetric-matrices/91843/5 "2022-12-19T15:49:31Z")

</div>

> [@corbat](#):
>
> (`product_3` is much slower and has a lot of allocations.)

Ah, sorry, I just noticed I didn’t properly replace one variable name and was accessing a global variable. `product_3` is actually fine regarding allocations.

Here is the speed comparison including `product_3`:

 ![btime_comparison_2](https://global.discourse-cdn.com/julialang/original/3X/3/0/30297d25f4f7f6ea278073798c9deede66f8c80c.png)

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [December 20, 2022, 4:24am UTC](https://discourse.julialang.org/t/faster-quadratic-expression-for-symmetric-matrices/91843/6 "2022-12-20T04:24:52Z")

</div>

I figured a symmetric quadratic form would be a fun topic for [a blog post](https://spmd.org/posts/optimizing_symmetric_quadratic_form/).

For a 200x200 test problem (an admittedly convenient size), I got the runtime for a symmetric matrix down to \<1 microsecond on my desktop.  
For comparison, `product_1` was 3.66 and 7 microseconds on dense and symmetric matrices, respectively, while `product_3` was 5.4 and 8.9.

`symv!` and a skew-symmetric version are a little more complicated because lining up computations is more important when you need to store them in a particular destination.

Perhaps Simon could follow this up with a write up on skew-symmetric matrix-vector multiplication (or perhaps I’ll choose that as my excuse to procrastinate someday); I’d be happy to accept a PR at [spmd\_blog](https://github.com/JuliaSIMD/spmd_blog).

---

<div class="post-metadata">

**Author:** ![corbat](https://avatars.discourse-cdn.com/v4/letter/c/e0b2c6/32.png) [@corbat](https://discourse.julialang.org/u/corbat)\
**Post date:** [December 20, 2022, 12:59pm UTC](https://discourse.julialang.org/t/faster-quadratic-expression-for-symmetric-matrices/91843/7 "2022-12-20T12:59:13Z")

</div>

Wow, thanks a lot for this very nice and detailed discussion of the problem.  
Even though your code is way above my skill level, the blog entry was a very interesting read for me.  
Looking forward to at some point use `@turbo` and letting “LoopModels” do all the magic 😉

---

<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:** [December 23, 2022, 6:04pm UTC](https://discourse.julialang.org/t/faster-quadratic-expression-for-symmetric-matrices/91843/8 "2022-12-23T18:04:20Z")

</div>

@Elrod ,As always, this is really great, 2 questions:

1. To make the comparison wit no allocations on both side, is there a way to pre allocated combination of `dsymv` and `matvec` product in `MKL`?
2. What about the case of large matrices with multi threaded code?

By the way, in case of using the same matrix repeatedly, what I do is as following:

1. Decompose the matrix.
2. Use the decomposition in the form of a matrix vector with `mul!()` with pre allocated buffer.
3. Use `dot()` on the result (Maybe `sum(abs2(), ...)` is better, anyway, there should be more efficient way).

---

<div class="post-metadata">

**Author:** ![blackeneth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/blackeneth/32/10353_2.png) [@blackeneth](https://discourse.julialang.org/u/blackeneth)\
**Post date:** [December 23, 2022, 6:49pm UTC](https://discourse.julialang.org/t/faster-quadratic-expression-for-symmetric-matrices/91843/9 "2022-12-23T18:49:39Z")

</div>

Great blog post. A few comments and questions:

**Typo?**

> CPI is cycles per instruction, the higher the value, the more clock cycles were required per instruction. We see that the execution stall rate was higher for product\_4.

I think you mean product\_3 here. You don’t define product\_4 until later.

**Question**

- What is a “rectangular loop”?

- could you setup an RSS feed for your blog?

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [December 23, 2022, 6:54pm UTC](https://discourse.julialang.org/t/faster-quadratic-expression-for-symmetric-matrices/91843/10 "2022-12-23T18:54:36Z")

</div>

a rectangular loop is when yo have nested loops where the bounds don’t interact with each other. (i.e. the shape the loop variables form is a hyper-rectangle)

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [December 27, 2022, 3:56am UTC](https://discourse.julialang.org/t/faster-quadratic-expression-for-symmetric-matrices/91843/11 "2022-12-27T03:56:58Z")

</div>

> [@blackeneth](#):
>
> I think you mean product\_3 here. You don’t define product\_4 until later.

Good catch, thanks.

> [@blackeneth](#):
>
> could you setup an RSS feed for your blog?

I’ve not used rss before, but I defined the global variables mentioned [here](https://franklinjl.org/syntax/rss/#global_configuration) in [this pr](https://github.com/JuliaSIMD/spmd_blog/commit/9378ad11670de8a5742cc4dfd5f4b55eb298c0c5).

Oscar is correct about a rectangular loop.  
The new LoopModels is going to support arbitrary polyhedra, and not just rectangles.  
That still doesn’t cover everything, e.g. many operations with sparse or ragged arrays aren’t Polyhedra.

Same with loops calling functions that have iterative algorithms, like `gcd`, inside of them. I plan on adding some support for those eventually, too…  
(Note that `gcd` currently does work with LoopVectorization, because it treats it as any old function, but if you manually inlined it, `while` loop and all, it won’t work)

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [December 27, 2022, 4:14am UTC](https://discourse.julialang.org/t/faster-quadratic-expression-for-symmetric-matrices/91843/12 "2022-12-27T04:14:18Z")

</div>

> [@RoyiAvital](#):
>
> - To make the comparison wit no allocations on both side, is there a way to pre allocated combination of `dsymv` and `matvec` product in `MKL`?

Sure:

```julia
julia> @btime product_8($x, $B)
  1.021 μs (0 allocations: 0 bytes)
560649.0855276876

julia> @btime product_2!($y, $x, $B)
  1.716 μs (0 allocations: 0 bytes)
560649.0855276876

julia> @btime product_2!($y, $x, $A)
  2.373 μs (0 allocations: 0 bytes)
560649.0855276876

```

Didn’t make much difference to minimum time.

> [@RoyiAvital](#):
>
> What about the case of large matrices with multi threaded code?

That should just be memory bound, but multithreading could certainly help by increasing the effective size of your local caches (and increasing the total memory bandwidth you have access to).  
Writing multithreaded versions would be a bit more cumbersome.  
At some point, I’ll improve Polyester’s local variable support, to give it an option to work more like `LoopVectorization`’s (which supports non-allocating reductions).  
Then we could `Polyester.@batch` the outer loop.

> [@RoyiAvital](#):
>
> Decompose the matrix.

I guess the point of this is that you can `dot(y,y)` or `sum(abs2, y)` instead of dotting two different vectors?  
`dot` is mostly limited by the rate at which you can load from the vectors…  
But I’m not sure how much this will help, especially if you can fuse the operations like I did.

But `x' * inv(S) * x` is also a really common operation, and that’s obviously a really good idea there.

Also, if you have the same matrix repeatedly, any chance you can multiply an entire matrix, instead of one vector at a time?  
`symm!` or `trmm!` are going to be way faster than a bunch of `symv!` or `trmv!`s.

---

<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:** [December 27, 2022, 6:04am UTC](https://discourse.julialang.org/t/faster-quadratic-expression-for-symmetric-matrices/91843/13 "2022-12-27T06:04:16Z")

</div>

@Elrod , Could you share the code for the non allocating versions?

My assumption is that in a mutli threaded scenario there will be gains in real world cases, hence it is a better representation. Maybe then different winner will emerge?

Indeed, when the one can, packing multiple mat vec operations into mat mat will reduce overhead. But sometimes the different vectors are available at different times, hence it is not viable.