# Hermitian version of three-argument dot

**URL:** <https://discourse.julialang.org/t/hermitian-version-of-three-argument-dot/106122>\
**Category:** General Usage\
**Tags:** linearalgebra, complex-numbers\
**Created:** [November 12, 2023, 5:39pm UTC](https://discourse.julialang.org/t/hermitian-version-of-three-argument-dot/106122 "2023-11-12T17:39:52Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![fph](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fph/32/17159_2.png) [@fph](https://discourse.julialang.org/u/fph)\
**Post date:** [November 12, 2023, 5:39pm UTC](https://discourse.julialang.org/t/hermitian-version-of-three-argument-dot/106122/1 "2023-11-12T17:39:53Z")

</div>

Is there a function to compute directly v^\* A v for a vector v and a matrix A?

Ideally, it should have the upsides of three-argument `dot`, like not allocating Av, and take advantage of symmetry, so when A is `Symmetric` or `Hermitian` and v is complex it should return a real value rather than a complex one.

```julia
v = randn(ComplexF64, 2, 1)
A = Symmetric(randn(2,2))
isreal(dot(v, A, v)) # often returns false

```

---

<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:** [November 12, 2023, 6:20pm UTC](https://discourse.julialang.org/t/hermitian-version-of-three-argument-dot/106122/2 "2023-11-12T18:20:14Z")

</div>

Not that I know of, besides `real(dot(v, A, v))` of course.

One option is to have `dot(x, A, y)` take a special branch if `x === y` (i.e. if they are the same object, which can be checked quickly), although of course it would still need to return a complex number (with zero imaginary part) for type stability when the vectors are complex. Since this is just an optimization, it could be added without much fanfare.

Another option is to have `norm(x, A)` return the `A`-norm \Vert x \Vert\_A = \sqrt{x^\* A x} for hermitian A (assumed positive-definite — it could throw a `DomainError` if A was indefinite whenever x^\* A x \< 0, since it’s too expensive to check positive-definiteness whenever you compute the norm).

---

<div class="post-metadata">

**Author:** ![fph](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fph/32/17159_2.png) [@fph](https://discourse.julialang.org/u/fph)\
**Post date:** [November 12, 2023, 6:30pm UTC](https://discourse.julialang.org/t/hermitian-version-of-three-argument-dot/106122/3 "2023-11-12T18:30:42Z")

</div>

Thanks, that makes sense.

Another option would be adding `norm2` to compute the square norm (we already have `abs2`, and `norm2` could be implemented in a more performant way than `norm` since one doesn’t need to check for overflow). Then `x'*A*x` could simply be a method `norm2(x, A)`. Has this function been considered before?

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [November 12, 2023, 6:39pm UTC](https://discourse.julialang.org/t/hermitian-version-of-three-argument-dot/106122/4 "2023-11-12T18:39:52Z")

</div>

Another (but slower) option would be to get more pedantic:

```julia
using IntervalArithmetic

intv = Interval.(v)

intd = dot(intv,A,intv)

contains_real(v) = contains_zero(imag(v))

```

Now,

```julia
julia> contains_real(intd)
true

julia> intd
[-0.414165, -0.414164] + [-8.32668e-17, 1.11023e-16]im

```

`contains_real` reliably returns `true`. Rounding towards a real can also be managed.

---

<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:** [November 12, 2023, 8:56pm UTC](https://discourse.julialang.org/t/hermitian-version-of-three-argument-dot/106122/5 "2023-11-12T20:56:13Z")

</div>

> [@fph](#):
>
> Another option would be adding `norm2` to compute the square norm (we already have `abs2`, and `norm2` could be implemented in a more performant way than `norm` since one doesn’t need to check for overflow)

Yes, though for an ordinary L2 norm squared you can simply do `sum(abs2, x)`.

---

<div class="post-metadata">

**Author:** ![dpo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpo/32/3335_2.png) [@dpo](https://discourse.julialang.org/u/dpo)\
**Post date:** [November 12, 2023, 9:14pm UTC](https://discourse.julialang.org/t/hermitian-version-of-three-argument-dot/106122/6 "2023-11-12T21:14:08Z")

</div>

As far as I understand, `sum(abs2, x)` doesn’t scale `x` prior to summing its components in order to avoid overflow in the way `norm()` does. Am I wrong here?

Also I would say that the computation of `x’ * A * x` is relevant also for symmetric/Hermitian indefinite `A`, not just for positive definite `A`.

---

<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:** [November 12, 2023, 9:21pm UTC](https://discourse.julialang.org/t/hermitian-version-of-three-argument-dot/106122/7 "2023-11-12T21:21:35Z")

</div>

> [@dpo](#):
>
> As far as I understand, `sum(abs2, x)` doesn’t scale `x` prior to summing its components in order to avoid overflow in the way `norm()` does. Am I wrong here?

Right, `sum(abs2, x)` does no scaling. Nor is it useful to do any scaling when you are computing the _squared_ norm.

It’s only if you are computing the non-squared norm, i.e. \sqrt{\sum |x\_k|^2}, that you want to use `norm(x)` instead of `sqrt(sum(abs2, x))` to avoid spurious over/underflow.

---

<div class="post-metadata">

**Author:** ![dpo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpo/32/3335_2.png) [@dpo](https://discourse.julialang.org/u/dpo)\
**Post date:** [November 12, 2023, 9:39pm UTC](https://discourse.julialang.org/t/hermitian-version-of-three-argument-dot/106122/8 "2023-11-12T21:39:19Z")

</div>

True.

---

<div class="post-metadata">

**Author:** ![mikmoore](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mikmoore/32/31109_2.png) [@mikmoore](https://discourse.julialang.org/u/mikmoore)\
**Post date:** [November 13, 2023, 4:24pm UTC](https://discourse.julialang.org/t/hermitian-version-of-three-argument-dot/106122/9 "2023-11-13T16:24:51Z")

</div>

Is there a algorithmic advantage to using `Hermitian`/`Symmetric` structure in `dot(v,A,v)`? It seems like one still needs N^2+N mul-adds regardless of the symmetry of `A`. Or is there an optimization I’m missing here?

_EDIT: the below post has shown that there is a 2x algorithmic advantage to some symmetry structures._

~~If you had some nice factorization of `A`, like `cholesky`, then you could save roughly half the computations. But this savings isn’t worth the cost of `cholesky` up-front. In the complex case, I suppose you could avoid computing the imaginary part of the second multiply to drop the cost from 4N^2+4N to 4N^2 + 2N real flops, but that would seem to be a marginal savings even at small sizes (24 vs 20 real flops at N=2 and more negligible at larger).~~

If what you care about is that the result is real because (structurally) it should be, consider just calling `real(dot(v,A,v))`.

---

<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:** [November 13, 2023, 5:58pm UTC](https://discourse.julialang.org/t/hermitian-version-of-three-argument-dot/106122/10 "2023-11-13T17:58:23Z")

</div>

> [@mikmoore](#):
>
> Is there an algorithmic advantage to using `Hermitian`/`Symmetric` structure in `dot(v,A,v)`?

Yes, you theoretically save about a factor of 2 because you only need to look at the upper triangle of A.

Let’s suppose that all the numbers are complex and you exploit no structure. Our `dot(v,A,v)` algorithm essentially computes the most straightforward:

\sum\_{i,j} \overline{v\_i} A\_{ij} v\_j = \sum\_i \overline{v\_i} \left( \sum\_j A\_{ij} v\_j \right)

which requires N^2 + N complex multiplications (6 flops each) and N(N - 1) + (N-1) = N^2 - 1 complex additions (2 flops each), for \approx 8N^2 flops. You can save O(N) flops by some tricks for `v === v`, e.g. special handling for the diagonal elements where you have a complex–real multiplication and a norm instead of two general complex–complex multiplies, but the leading term is still 8N^2. (This is just a naive flop count, of course, not counting FMA or SIMD.)

For Hermitian A, however, things simplify to

\sum\_{i=1}^N |v\_i|^2 A\_{ii} + 2 \sum\_{i=1}^N \operatorname{Re}\left(\overline{v\_i} \left[\sum\_{j=i+1}^N A\_{ij} v\_j \right]\right)

which costs N `abs2` computations (3 flops each) + N real multiplies + N-1 real additions (1 flop each) for the first sum, along with N(N-1)/2 general complex multiplies (6 flops each) + N complex multiplies where you only compute the real part (3 flops each) + N(N-1)/2 - N complex additions (2 flops each) + N-1 real additions (1 flop each) for the second sum, + 1 real add and 1 real multiply (2 flops) to put it all together. The leading term is 3N^2 + N^2 = 4N^2 flops, **asymptotically saving a factor of 2.**

If all of the numbers are real, then the savings in flops are still asymptotically a factor of 2 (the leading term goes from \approx 2N^2 flops to \approx N^2 flops). Of course, for large matrices these BLAS2-type operations are typically memory-bound and not compute-bound, but there should still be memory-bandwidth savings from only looking at half of A.

---

<div class="post-metadata">

**Author:** ![mikmoore](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mikmoore/32/31109_2.png) [@mikmoore](https://discourse.julialang.org/u/mikmoore)\
**Post date:** [November 13, 2023, 6:44pm UTC](https://discourse.julialang.org/t/hermitian-version-of-three-argument-dot/106122/11 "2023-11-13T18:44:05Z")

</div>

Thanks. So the short version is that for Hermitian `A` you can compute something algorithmically like `2*real(dot(v,triu(A,1),v)) + dot(v,Diagonal(A),v)`, which does indeed induce the sort of sparsity one can exploit to shrink the leading coefficient on the flop count.

These exploits appear to require that `x == y` (in `dot(x,A,y)`), so it’s quite a pun to call this `dot` any more. As someone pointed out above, it’d better be described as a mode of `norm`.

_EDIT: as pointed out below, this is merely a performance enhancement so isn’t “wrong” to include it in `dot` when we can observe that the vector arguments are `===`. For type stability, the result would still be complex-typed for complex arguments but we could ensure the imaginary part is zero. Although `===` is cheap for `Array`, it would have more significant consequences if this method ended up applying to bitstypes like `StaticArrays.SVector`._

---

<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:** [November 13, 2023, 6:51pm UTC](https://discourse.julialang.org/t/hermitian-version-of-three-argument-dot/106122/12 "2023-11-13T18:51:16Z")

</div>

> [@mikmoore](#):
>
> These exploits appear to require that `x == y` (in `dot(x,A,y)`), so it’s quite a pun to call this `dot` any more. As someone pointed out above, it’d better be described as a mode of `norm`.

Right, but that being said there’s nothing wrong with writing a method of `dot(x, A::Hermitian, y)` that checks `x === y` (not `x == y`, which is slower) and goes to an optimized implementation for that case.