# A \* B \* A' is not hermitian, even when B is

**URL:** <https://discourse.julialang.org/t/a-b-a-is-not-hermitian-even-when-b-is/70611>\
**Category:** General Usage\
**Tags:** linearalgebra\
**Created:** [October 29, 2021, 11:08am UTC](https://discourse.julialang.org/t/a-b-a-is-not-hermitian-even-when-b-is/70611 "2021-10-29T11:08:48Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![jlbosse](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jlbosse/32/11274_2.png) [@jlbosse](https://discourse.julialang.org/u/jlbosse)\
**Post date:** [October 29, 2021, 11:08am UTC](https://discourse.julialang.org/t/a-b-a-is-not-hermitian-even-when-b-is/70611/1 "2021-10-29T11:08:48Z")

</div>

The title essentially says it all: If I create a hermitian matrix `A` with

```julia
julia> R = randn(10,10)
julia> A = R * R' # definitely hermitian and positive semidefinite

```

and another arbitrary matrix `B`

```julia
julia> B = rand(7, 10)

```

then

```julia
julia> ishermitian(B * A * B')
false

```

even though, mathematically, it should be hermitian. I suspect the problem is simply some issues with floating point arithmetic. Is there a fast method that guarantees that the result is hermitian?

---

<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:** [October 29, 2021, 11:14am UTC](https://discourse.julialang.org/t/a-b-a-is-not-hermitian-even-when-b-is/70611/2 "2021-10-29T11:14:11Z")

</div>

> [@jlbosse](#):
>
> I suspect the problem is simply some issues with floating point arithmetic. Is there a fast method that guarantees that the result is hermitian?

Yes, it’s just roundoff errors. You can use `Hermitian(B * A * B')`.

---

<div class="post-metadata">

**Author:** ![mcabbott](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mcabbott/32/6603_2.png) [@mcabbott](https://discourse.julialang.org/u/mcabbott)\
**Post date:** [October 29, 2021, 12:04pm UTC](https://discourse.julialang.org/t/a-b-a-is-not-hermitian-even-when-b-is/70611/3 "2021-10-29T12:04:43Z")

</div>

Since `B * A * B'` is initially one function call to `*`, it’s possible that if `A::Hermitian`, this could propagate? I think it needs to check that both sides are the same `B` though, so it can’t be type-stable, and is probably not a great idea:

```julia
julia> A = Hermitian(R * R');

julia> @less B * A * B' # 3-matrix method, calls _tri_matmul(A,B,C) ... on Julia >= 1.7?

julia> Base.:*(A::AbstractMatrix, B::Hermitian, C::Adjoint{<:Any, <:AbstractMatrix}) = A === parent(C) ? Hermitian(A * (B * C)) : LinearAlgebra._tri_matmul(A,B,C)

julia> B * A * B'
7×7 Hermitian{Float64, Matrix{Float64}}:
 15.6641 14.2528 22.1946 8.88175 19.089 12.2722 14.4311
 14.2528 29.7732 33.362 14.4554 22.016 19.0204 19.188
...

julia> @code_warntype B * A * B'
...
Body::Union{Hermitian{Float64, Matrix{Float64}}, Matrix{Float64}}

```

---

<div class="post-metadata">

**Author:** ![mtfishman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mtfishman/32/30755_2.png) [@mtfishman](https://discourse.julialang.org/u/mtfishman)\
**Post date:** [October 29, 2021, 1:52pm UTC](https://discourse.julialang.org/t/a-b-a-is-not-hermitian-even-when-b-is/70611/4 "2021-10-29T13:52:28Z")

</div>

I’m sure this has been brought up in the past, but shouldn’t `ishermitian` have optional tolerance cutoffs?

Internally in code you might want to check if a matrix is approximately Hermitian and if so wrap it in `Hermitian` to call specialized algorithms. Though I guess you can do it yourself with `isapprox(M, M'; ...)`.

---

<div class="post-metadata">

**Author:** ![jishnub](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jishnub/32/33620_2.png) [@jishnub](https://discourse.julialang.org/u/jishnub)\
**Post date:** [October 29, 2021, 2:03pm UTC](https://discourse.julialang.org/t/a-b-a-is-not-hermitian-even-when-b-is/70611/5 "2021-10-29T14:03:19Z")

</div>

Was brought up 10 days back: [https://github.com/JuliaLang/julia/pull/42707](https://github.com/JuliaLang/julia/pull/42707)

---

<div class="post-metadata">

**Author:** ![jlbosse](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jlbosse/32/11274_2.png) [@jlbosse](https://discourse.julialang.org/u/jlbosse)\
**Post date:** [October 29, 2021, 5:06pm UTC](https://discourse.julialang.org/t/a-b-a-is-not-hermitian-even-when-b-is/70611/6 "2021-10-29T17:06:52Z")

</div>

I suppose just wrapping it in `Hermitian(B * A * B')` is easiest for now. If I get around to it (and need the extra performance) it might be advantageous to write an extra function that does this tri-product `hermitiansandwich(A, B) = ...` that has hermiticity baked in.

---

<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:** [October 29, 2021, 5:08pm UTC](https://discourse.julialang.org/t/a-b-a-is-not-hermitian-even-when-b-is/70611/7 "2021-10-29T17:08:08Z")

</div>

> [@mtfishman](#):
>
> Internally in code you might want to check if a matrix is approximately Hermitian and if so wrap it in `Hermitian` to call specialized algorithms

In most cases, in your own code you should _know_ whether a matrix is supposed to be Hermitian in exact arithmetic. If users are passing you arbitrary matrices, it’s reasonable to ask them to indicate this by using the `Hermitian` type if possible.

The challenge with using `isapprox(A, A', rtol=???)` to decide whether you can use some Hermitian-specialized algorithm is that you then need to do some careful error analysis to decide what tolerance to use.

---

<div class="post-metadata">

**Author:** ![mtfishman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mtfishman/32/30755_2.png) [@mtfishman](https://discourse.julialang.org/u/mtfishman)\
**Post date:** [October 29, 2021, 6:18pm UTC](https://discourse.julialang.org/t/a-b-a-is-not-hermitian-even-when-b-is/70611/8 "2021-10-29T18:18:23Z")

</div>

Yes, I suppose that’s true. We have use cases in tensor networks where we reduce a complicated set of operations (for example, a series of tensor contractions, eigendecomositions, etc.) to a matrix which is supposed to be Hermitian but may not be in practice due to numerical errors or coding errors. So it generalizes `B * A * B'` to cases where it is much less obvious that the resulting matrix ends up Hermitian. But again, these cases can be checked manually with `isapprox` and doesn’t need to be inside `ishermitian`.

---

<div class="post-metadata">

**Author:** ![jlapeyre](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jlapeyre/32/4514_2.png) [@jlapeyre](https://discourse.julialang.org/u/jlapeyre)\
**Post date:** [November 1, 2021, 1:17pm UTC](https://discourse.julialang.org/t/a-b-a-is-not-hermitian-even-when-b-is/70611/9 "2021-11-01T13:17:11Z")

</div>

I wrote a package a while back to try to treat these questions in a systematic way: [IsApprox.jl](https://github.com/jlapeyre/IsApprox.jl)

```julia
julia> using IsApprox

julia> R = randn(10,10); A = R * R';

julia> B = rand(7, 10);

julia> M = B * A * B';

julia> IsApprox.ishermitian(M)
false

julia> IsApprox.ishermitian(M, Equal())
false

julia> IsApprox.ishermitian(M, Approx())
true

```