# Faster way for matrix multiplication including a Hermitian matrix with special properties

**URL:** <https://discourse.julialang.org/t/faster-way-for-matrix-multiplication-including-a-hermitian-matrix-with-special-properties/103177>\
**Category:** New to Julia\
**Tags:** quantum, linearalgebra, optimization\
**Created:** [August 25, 2023, 7:31am UTC](https://discourse.julialang.org/t/faster-way-for-matrix-multiplication-including-a-hermitian-matrix-with-special-properties/103177 "2023-08-25T07:31:37Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![ahmed](https://avatars.discourse-cdn.com/v4/letter/a/8c91f0/32.png) [@ahmed](https://discourse.julialang.org/u/ahmed)\
**Post date:** [August 25, 2023, 7:31am UTC](https://discourse.julialang.org/t/faster-way-for-matrix-multiplication-including-a-hermitian-matrix-with-special-properties/103177/1 "2023-08-25T07:31:37Z")

</div>

I have a simple question that I hope can make certain kinds of matrix multiplications faster in my code.  
Is there a way to make use of matrix multiplications that are of the form `phi*P*adjoint(phi)` to make the computation faster? where P is a low rank Hermitian matrix, specifically it has eigen values of only ones and zeros (half of the eigenvalues are ones). phi is in general a complex matrix (not a square matrix) but its columns are orthogonal.  
I would think that the combination of the Hermiticity of P, its low rank and the orthogonality of the columns vectors of phi can work to optimize these kinds of products. It seems like the standard linear algebra package doesn’t make a use of that as illustrated in this example

```julia
using LinearAlgebra, BenchmarkTools
A1 = rand(ComplexF64, 10, 10);
A = A1 + adjoint(A1); #make a hermitian matrix
phi = rand(ComplexF64, 1000, 10);
function test(phi, A) # test the performance when A is Hermitian 
	return phi*A*adjoint(phi)
end 
@btime test(phi, A)

function test2(phi, A1) # test the performance where A1 is not Hermitian 
	return phi*A1*adjoint(phi)
end
@btime test2(phi, A1)

```

The time of the calculations with a Hermitian matrix and a not Hermitian matrix comes out to be the same where if The symmetry is used; the Hermitian matrix calculation should be roughly half the time.  
Any insight about this problem would be appreciated, Thanks!

---

<div class="post-metadata">

**Author:** ![jd-foster](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jd-foster/32/35824_2.png) [@jd-foster](https://discourse.julialang.org/u/jd-foster)\
**Post date:** [August 25, 2023, 10:43am UTC](https://discourse.julialang.org/t/faster-way-for-matrix-multiplication-including-a-hermitian-matrix-with-special-properties/103177/2 "2023-08-25T10:43:50Z")

</div>

It could be worthwhile explicitly using the `Hermitian` constructor: see [Linear Algebra · The Julia Language](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#Special-matrices)

---

<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:** [August 25, 2023, 3:47pm UTC](https://discourse.julialang.org/t/faster-way-for-matrix-multiplication-including-a-hermitian-matrix-with-special-properties/103177/3 "2023-08-25T15:47:56Z")

</div>

As the previous person suggested, you should wrap the matrix in `Hermitian(A)` if you want Julia to know it’s Hermitian. There are only a handful of functions where Julia speculatively checks for properties like Hermitian or triangular or diagonal to decide what algorithm to use. That said, I don’t think it will help here in any case.

If your `A` has a symmetric factorization (such as a Cholesky or LDL^\top factorization), then you can get the first of your two multiplies for half cost and change your overall computation (e.g., to `(phi * L) * D * (phi * L)'`, for the LDL^\top case). But unless `A` is small relative to `phi` (like in your example), this factorization may be more expensive than the savings it offers. In any case, it’s something you should apply manually.

So, of the two multiplies, you can save half on the first if `A` is symmetrically factorized. Further, you could\* save half on the second because the result is Hermitian.

\*_HOWEVER_, I am [not aware of a routine in BLAS](https://stackoverflow.com/questions/74545160/how-to-tell-blas-that-output-matrix-is-symmetric) (or other libraries or Julia) for a matrix multiply function that only computes half of the output matrix. This would be necessary for the savings on the second multiply. So unless someone else does, those savings are off the table. This leaves you with either no savings (general `A`) or 25% total savings (suitably factorized `A`).

That’s my assessment, anyway. Maybe someone else knows of a suitable half-output matrix multiply? Or some radically different approach?

---

<div class="post-metadata">

**Author:** ![ahmed](https://avatars.discourse-cdn.com/v4/letter/a/8c91f0/32.png) [@ahmed](https://discourse.julialang.org/u/ahmed)\
**Post date:** [August 25, 2023, 6:29pm UTC](https://discourse.julialang.org/t/faster-way-for-matrix-multiplication-including-a-hermitian-matrix-with-special-properties/103177/4 "2023-08-25T18:29:53Z")

</div>

Thanks @jd-foster and @mikmoore for your replies. Just to clarify if ishermitian(A) returns true, is it also still required to do Hermitian(A) for the linearalgebra package to specialize that A is Hermitian? (I’d like to know this for future reference even if it doesn’t help me in this specific example). It seems like by itself the LinearAlgebra package doesn’t speed the calculation even if I explicitly add Hermitian(A).  
@mikmoore wouldn’t what you describe also work if A has an eigenvalue decomposition?

---

<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:** [August 25, 2023, 7:07pm UTC](https://discourse.julialang.org/t/faster-way-for-matrix-multiplication-including-a-hermitian-matrix-with-special-properties/103177/5 "2023-08-25T19:07:42Z")

</div>

For some functions, Julia will check the structure up-front to decide what algorithm to use. `edit` (or the macro version `@edit`) is useful if you want to answer this question for some specific function. See `edit(\, (AbstractMatrix, AbstractVector))` and `edit(eigen, (AbstractMatrix,))` for examples of functions that do.

But I wouldn’t count on this in general. It’s best to wrap with `Hermitian` manually if you know your matrix has that structure and you think it might have a specialized algorithm.

> [@ahmed](#):
>
> wouldn’t what you describe also work if A has an eigenvalue decomposition?

Not really. It’s a Hermitian matrix, so the eigendecomposition is `V*D*V'` for a unitary `V` and diagonal `D`. So we could implement your product as

```julia
phi_V = phi * V
result = phi_V * D * phi_V'

```

Let’s use your example sizes `size(phi) == (1000,10)` and `size(A) == size(V) == size(D) == (10,10)`. `V` likely doesn’t have a sparsity pattern we could exploit to multiply more efficiently, so `phi*V` requires `1000*10*10` multiply-adds to produce a `size == (1000,10)` result. Then `phi_V * D` requires `1000*10` multiplies because `D` is diagonal. Then the final product requires `1000*10*10` multiply-adds, which we could reduce to roughly `1000*10*(10+1)/2` _if we had a function to only compute half the result_ (which I don’t believe we do, nor has BLAS felt the need to have one).

Notably, `phi*V` is equally expensive as `phi*A`. The fact that the eigendecomposition method requires the application of the diagonal matrix means it actually requires _slightly more_ work than `phi * A * phi'`.

The reason that a `L*L'` or `L*D*L'` factorization helps is that `L` is triangular. Multiplying `phi*L` needs only half the computation of `phi*V` or `phi*A`. But if you don’t already have `A` appropriately factorized, factorizing it may cost more than it saves.

* * *

The short summary is that you could theoretically shave 25% of the calculation with a suitably factorized `A` and another 25% if you can find an optimized half-output matrix multiply (if you write it yourself and aren’t a BLAS expert, it will likely be _much_ slower than producing the full output with an optimized version). Go ahead and do what you can with what you have handy, but I expect it isn’t worth getting too tangled up in trying to optimize this.

---

<div class="post-metadata">

**Author:** ![ahmed](https://avatars.discourse-cdn.com/v4/letter/a/8c91f0/32.png) [@ahmed](https://discourse.julialang.org/u/ahmed)\
**Post date:** [August 25, 2023, 10:39pm UTC](https://discourse.julialang.org/t/faster-way-for-matrix-multiplication-including-a-hermitian-matrix-with-special-properties/103177/6 "2023-08-25T22:39:20Z")

</div>

Thanks @mikmoore for your explanation. I’m not sure that my A has a symmetric factorization. (edit: an eigenvalue decomposition of A will give a diagonal matrix of ones and zeros. So I think this means that it doesn’t have a symmetric factorization as it is not positive definite) This means, however, that it is low rank, and that also should factor in speeding up the multiplication. I have edited the main post to include this information.  
It’s kind of a bummer that there is no built in optimization for these kind of products 🙂 especially I’d think they’re pretty common. And yes, I’m not really the qualified person that can take on this task of developing a package for that 🙂

---

<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:** [August 26, 2023, 2:38pm UTC](https://discourse.julialang.org/t/faster-way-for-matrix-multiplication-including-a-hermitian-matrix-with-special-properties/103177/7 "2023-08-26T14:38:40Z")

</div>

We do have a specialization for `dot(x,A,y) == x'*A*y` when `x` and `y` are vectors, but that doesn’t apply to your case.

If you eigendecomposition has zeros, your matrix is rank deficient. Since the matrix is Hermitian and the nonzero eigenvalues are positive, it is positive semidefinite. This means that you can write your `A` as `B'*B` where `size(B) == (R,N)` with `R == rank(A) < N`. This will reduce the cost to `R/N` of the naive cost on both multiplies.

Finding a low rank factorization is something you can do from an eigendecomposition. Just drop the (near-)zero terms. But you might already have `A` in a low rank form based on how you computed it in the first place.

Also your final product `phi * A * phi'` is very low rank. If you need all the entries, then fine. But if you’re using it in further calculations, keeping _it_ low-rank factorized is likely much more efficient.