I'm trying to simulate a problem in physics for which I require a Unitary operator in a Hilbert space with an inner product defined as transpose(x*)x. Given two orthogonal column vectors, I want to generate more orthonormal vectors. Here is the way I tried approaching this problem. I randomly generate complex vectors and subtract their projections onto other already available orthogonal vectors Similar to this. And then I check the norm with respect to the inner product (InProd). Here is an attempt to implement this using python.
def stinemod():
comped = [[1,0,0,0,0],[0,1,0,0,0]]
d = len(comped[0])
r = len(comped)
while(r<d):
randr = np.random.rand(d)
randc = np.random.rand(d)
vr = randr + 1j*randc
vo = vr
for v in comped:
vo = vo - (InProd(vr,v)*np.array(v)/InProd(v,v))
k = 1e-10
if(InProd(vo,vo)<k*InProd(vr,vr)):
pass
else:
r = r+1
comped.append(np.array(vr)/InProd(vr,vr))
return(np.transpose(comped))
But on running this code and checking unitarity using,
A = stinemod()
print(abs(np.matmul(np.transpose(np.conj(A)),A)))
Output:
[[1. 0. 0.28003392 0.24068132 0.1977418 ]
[0. 1. 0.53992755 0.24199218 0.06786818]
[0.28003392 0.53992755 0.58108559 0.29561698 0.23971144]
[0.24068132 0.24199218 0.29561698 0.21599542 0.18374313]
[0.1977418 0.06786818 0.23971144 0.18374313 0.21586778]]
I get an output suggesting that it is not Unitary which means columns are not orthonormal. I can't seem to figure out what the mistake in here is.