Unnecessary matrix products

While looking at !62 (merged) I saw that we are making more matrix products than necessary. I used the following code to reproduce (requires !62 (merged)):

# %%
from collections import Counter
from pymablock import block_diagonalize
from pymablock.series import AlgebraElement, zero
from matplotlib import pyplot as plt
import numpy as np

def solve_sylvester(A):
    return AlgebraElement(f"S({A})")

H_tilde, *_ = block_diagonalize(
    [[[AlgebraElement("(H_{0,AA}"), zero], [zero, AlgebraElement("H_{0,BB}")]],
     [[AlgebraElement("H_{1,AA}"), AlgebraElement("H_{1,AB}")],
      [AlgebraElement("H_{1,BA}"), AlgebraElement("H_{1,BB}")]]],
    solve_sylvester=solve_sylvester,)

# H_tilde, *_ = block_diagonalize(
#     [[[AlgebraElement("(H_{0,AA}"), zero], [zero, AlgebraElement("H_{0,BB}")]],
#      [[zero, AlgebraElement("H_{1,AB}")],
#       [AlgebraElement("H_{1,BA}"), zero]]],
#     solve_sylvester=solve_sylvester,)
# %%
AlgebraElement.log = []
multiplication_counts = []
for order in range(10):
    H_tilde[0, 0, order]
    multiplication_counts.append(Counter(call[1] for call in AlgebraElement.log)["__mul__"])
# %%
plt.plot(range(1, 10), np.diff(multiplication_counts))
plt.xlabel("Order")
plt.ylabel("Multiplication count")
# %%
AlgebraElement.log = []
# %%
H_tilde[0, 0, 2]
AlgebraElement.log

The 2nd order expression becomes

((((-(-S(H_{1,AB}))) \times H_{1,BA}) + (H_{1,AB} \times (-((-S(H_{1,AB}))^\dagger)))) + ((((-(-S(H_{1,AB}))) \times (H_{1,AB}^\dagger)) + (((-(-S(H_{1,AB}))) \times (H_{1,AB}^\dagger))^\dagger)) / -2))

This uses 4 matrix products, however direct inspection shows that all 4 are either identical (but don't use the Hermiticity of H' or Hermitian conjugates of each other