I wrote some code to implement the modified Gram Schmidt process. When I tested it on real matrices, it is correct. However, when I tested it on complex matrices, it went wrong.

I believe my code is correct by doing a step by step check. Therefore, I wonder if there are numerical reasons why the modified Gram Schmidt process fails on complex vectors.

Following is the code:

import numpy as np

def modifiedGramSchmidt(A):
    """
    Gives a orthonormal matrix, using modified Gram Schmidt Procedure
    :param A: a matrix of column vectors
    :return: a matrix of orthonormal column vectors
    """
    # assuming A is a square matrix
    dim = A.shape[0]
    Q = np.zeros(A.shape, dtype=A.dtype)
    for j in range(0, dim):
        q = A[:,j]
        for i in range(0, j):
            rij = np.vdot(q, Q[:,i])
            q = q - rij*Q[:,i]
        rjj = np.linalg.norm(q, ord=2)
        if np.isclose(rjj,0.0):
            raise ValueError("invalid input matrix")
        else:
            Q[:,j] = q/rjj
    return Q

Following is the test code:

import numpy as np

# If testing on random matrices:
# X = np.random.rand(dim,dim)*10  + np.random.rand(dim,dim)*5 *1j
# If testing on some good one
v1 = np.array([1, 0, 1j]).reshape((3,1))
v2 = np.array([-1, 1j, 1]).reshape((3,1))
v3 = np.array([0, -1, 1j+1]).reshape((3,1))
X = np.hstack([v1,v2,v3])
Y = modifiedGramSchmidt(X)
Y3 = np.linalg.qr(X, mode="complete")[0]
if np.isclose(Y3.conj().T.dot(Y3), np.eye(dim, dtype=complex)).all():
    print("The QR-complete gives orthonormal vectors")
if np.isclose(Y.conj().T.dot(Y), np.eye(dim, dtype=complex)).all():
    print("The Gram Schmidt process is tested against a random matrix")
else:
    print("But My modified GS goes wrong!")
    print(Y.conj().T.dot(Y))

Update

The problem is that I implemented a algorithm designed for inner product linear in first argument

whereas I thought it were linear in second argument.

Thanks @landogardner