How to Generate a Random Complex Symmetric Unitary Matrix in Python?

I am trying to optimize a complex symmetric unitary matrix 𝐴 with dimensions N*N using Python.

First of all, how to generate matrix 𝐴 that should satisfies the two conditions:

𝐴=𝐴^{T} (symmetric constraint), where 𝐴^{T} is the transpose of 𝐴

𝐴 . 𝐴^{H}= I (Unitary constraint), where 𝐴^{H} is the Hermitian transpose of 𝐴 and I is the identity matrix.

Secondly, I will use it to maximize a matrix multiplication of matrix B with N*N dimensions, so how to perform this matrix multiplication while 𝐴 remains both symmetric and unitary.

import random
import numpy as np
import cmath
import math

N = 4  # Matrices dimensions 

B=np.array([[0.09+0.76j, 0.67+0.08j, 0.28+0.58j],
            [0.47+0.82j, 0.85+0.07j, 0.13+0.49j],
            [0.12+0.8j , 0.35+0.31j, 0.76+0.85j]])                  

result =  A @ B

First find the unitary matrix U (the bulk of the stuff below). Then V=U^TU will be symmetric and unitary.

The harder part is getting the unitary matrix U. A matrix is unitary if and only if its columns are orthonormal. (Dot products of different columns are 0; norm of each column is 1.) So you can construct such a matrix by the Gram-Schmidt process.

Imagine each column of your matrix is a vector. Dealing with each column in turn you (a) subtract the projections of previous column vectors; then (b) normalise this column. As long as you have independent columns (i.e. your original matrix was invertible) then this will produce orthonormal columns and hence a unitary matrix.

The following will produce a unitary matrix from a given matrix A. A unitary matrix is one for which the conjugate of its transpose is its inverse (see https://en.wikipedia.org/wiki/Unitary_matrix ). It tests this at the end by multiplication to give the identity matrix (up to the usual floating-point accuracy).

You can improve this by (a) checking matrix A is invertible first (e.g. by calculating its determinant) and (b) formatting the output a bit more neatly.

import numpy as np

def makeUnitary( A ):
    N = A.shape[0]                          # number of columns
    U = A.copy()
    # Gram-Schmidt process
    U[:,0] /= np.linalg.norm( U[:,0] )      # normalise first column
    for i in range( 1, N ):                 # make successive columns orthogonal
        for j in range( i ):                #    by subtracting projections of previous columns
            U[:,i] -= U[:,i].dot( U[:,j].conj() ) * U[:,j]
        U[:,i] /= np.linalg.norm( U[:,i] )  # normalise column
    return U


A=np.array([[0.09+0.76j, 0.67+0.08j, 0.28+0.58j],
            [0.47+0.82j, 0.85+0.07j, 0.13+0.49j],
            [0.12+0.8j , 0.35+0.31j, 0.76+0.85j]])               
U = makeUnitary( A )                                  # a unitary matrix
V = U.T @ U                                           # symmetrised unitary matrix


print( "Original matrix:\n", A )
print( "\nSymmetric unitary matrix constructed from it:\n", V )
print( "\nCheck V.V(conj)T (should be identity matrix):\n", V @ V.conj().T )

Output:

Original matrix:
 [[0.09+0.76j 0.67+0.08j 0.28+0.58j]
 [0.47+0.82j 0.85+0.07j 0.13+0.49j]
 [0.12+0.8j  0.35+0.31j 0.76+0.85j]]

Symmetric unitary matrix constructed from it:
 [[-0.77181963+0.51542139j  0.15837321+0.05726809j -0.30709606+0.12635484j]
 [ 0.15837321+0.05726809j  0.40398494-0.53221372j -0.43217941-0.581725j  ]
 [-0.30709606+0.12635484j -0.43217941-0.581725j    0.55047024-0.24804427j]]

Check V.V(conj)T (should be identity matrix):
 [[1.00000000e+00-2.05296819e-18j 3.93579786e-16-6.38754713e-16j
  6.37756498e-17-6.66363810e-16j]
 [3.93579786e-16+6.24994617e-16j 1.00000000e+00-1.58327863e-18j
  1.09796762e-16-1.28773584e-16j]
 [6.37756498e-17+6.55406471e-16j 1.09796762e-16+1.49124945e-16j
  1.00000000e+00-8.89174271e-19j]]

A unitary matrix can be diagonalised (i.e. written as Q^-1.diag.Q). Powers of U can then be computed quickly simply by taking powers of the elements in the diagonal matrix.