eigh incompatible with multidirectional UTP?

It seems that eigh is inherently incompatible with multidirectional UTP in the presence of repeated eigenvalues. The Python implementation seems correct wrt. the mathematical theorems, but can only work if the UTP is unidirectional. Here is a case reproducing the problem:

from algopy import *

D = 3
P = 2
N = 6

A = UTPM(numpy.random.random((D, P, N, N)))

# Ensure all directions have the same constant value in A (and hence later Z1)
A.data[0,:,:,:] = A.data[0,0,:,:]

# Make an orthonormal basis Z1
[Z1, _] = qr(A)

# Orthogonal basis Z1 has same constant value in all directions
numpy.testing.assert_almost_equal((Z1.data[0,:,:,:]-Z1.data[0,0,:,:]), 0)

# Make fresh eigenvalues with some repetitions
W1 = zeros((N, N), dtype=A)
for d in range(D):
    for p in range(P):
        for i in range(N):
            W1.data[d, p, i, i] = (1 + (d*(p + 1) + round(i / 2))) * 100
W1.data[0,:,:,:] = W1.data[0,0,:,:]

# Matrix A1 is the one we want to factorize into eigenvalues and vectors
A1 = dot(Z1, dot(W1, Z1.T))

# Factorize A1 as w2 and Z2
[w2, Z2] = eigh(A1)

# Z2 is indeed orthogonal as required
numpy.testing.assert_almost_equal(dot(Z2, Z2.T) - numpy.eye(N), 0)

# Reconstruction A2 of A1 via its eigendecomposition succeeds - as it should!
W2 = zeros_like(A)
W2[:N, :N] = diag(w2)
A2 = dot(Z2, dot(W2, Z2.T))
numpy.testing.assert_almost_equal(A1.data-A2.data, 0)

# But Z2 does not have the same constant value in each direction!
print("Z2 does not have the same constant value in all directions:\n\n",
        Z2.data[0,:,:,:], "\n")
numpy.testing.assert_almost_equal((Z2.data[0,:,:,:]-Z2.data[0,0,:,:]), 0)

Could you please confirm whether this issue cannot be solved for multidirectional UTP, i.e. symmetric eigendecompositions containing repeated eigenvalues can only ever work in the unidirectional case?

b45ch1 / algopy

eigh incompatible with multidirectional UTP? #54