import numpy as np
import scipy as sp
import scipy.optimize as so

numMatrices = 1
n = 3

def hermitian(A, **kwargs):
    return np.transpose(A,**kwargs).conj()

H = hermitian

def matrixProduct(a, lowerTriag):
    # We specified everything below the diagonal
    # as known numbers to reduce the number of
    # unknowns to match the number of equations
    aMod = lowerTriag + np.triu(np.reshape(a,(n,n)))
    ret = np.einsum('ij,jk', H(aMod), aMod)-np.eye(np.shape(aMod)[0])
    return ret.flatten()

lowerTriang = []
finalMatrix = []

for mat in xrange(numMatrices):
    lowerTriang.append(np.tril(np.random.random((n,n)),-1))

    x0 = lowerTriang[-1] + np.triu(np.random.random((n,n)))

    finalMatrix.append(so.fsolve(matrixProduct,
                                 x0.flatten(),
                                 args=(lowerTriang[-1]),
                                 maxfev=100000))

for mat in finalMatrix:
    newMat = np.reshape(mat,(n,n))
    print np.dot(H(newMat),newMat)