Orthogonal, Mean-Invariant Matrix Generation Between Python and Matlab

Question

I needed a program that produces a random orthogonal matrix, which also preserves the mean and covariance of the columns of a matrix. In other words, the matrix should have the property of any orthogonal matrix A*A.T = I and also have the mean-preserving property of A*1n = 1n, where 1n is a column of the scalar 1 of size 1xn. The technique to do that is pretty simple:

1. Produce a random matrix
2. Concatenate a column of 1s as the first column
3. Perform Gram-Schmidt orthonormalization on the matrix
4. Do steps 1, 2 and 3 for another matrix
5. Multiply these two matrices

I have a program in Matlab that sucessfully does that:

function M = GenerateROM(SeedValue, p)
rng(SeedValue);
Z1 = [ones(p, 1) randn(p, p-1)];
M1 = GramSchmidtOrthonorm(Z1);
rng(SeedValue+2);
Z2 = [ones(p, 1) randn(p, p-1)];
M2 = GramSchmidtOrthonorm(Z2);
M = M1 * M2';

function M = GramSchmidtOrthonorm(Z)
[p, col] = size(Z);
Y = [];
M = [];
for i = 1:p
    v = Z(:, i);
    u = v;
    for j = 1:(i-1)
        y = Y(:, j);
        u = u - (y' * v) / (y' *y) *y;
    end
    Y = [Y u];
    M = [M u/sqrt(u' *u)];
end

The above code does indeed satisfy the two properties I need.

I want to write similar code in Python. I tried the following:

import numpy as np
from scipy.stats import ortho_group

def gram_schmidt(A):
    """Orthogonalize a set of vectors stored as the columns of matrix A."""
    # Get the number of vectors.
    n = A.shape[1]
    for j in range(n):
        # To orthogonalize the vector in column j with respect to the
        # previous vectors, subtract from it its projection onto
        # each of the previous vectors.
        for k in range(j):
            A[:, j] -= np.dot(A[:, k], A[:, j]) * A[:, k]
        A[:, j] = A[:, j] / np.linalg.norm(A[:, j])
    return A

onesVector = np.ones((aRows,1))

random1 = np.random.normal(size=(aRows,aRows-1))
Z1 = np.hstack((onesVector, random1))
M1 = gram_schmidt(Z1)

random2 = np.random.normal(size=(aRows,aRows-1))
Z2 = np.hstack((onesVector, random2))
M2 = gram_schmidt(Z2)

A = M1@M2
print('This is A*A^T: ', np.around([email protected]))
print('This is A*1: ', np.around(A@onesVector))

Although this Python code satisfies the 1st property and gives me the identity matrix, it does not satisfy the 2nd property as shown below:

This is A*A^T:  [[1.00 -0.00 0.00 ... -0.00 0.00 -0.00]
[-0.00 1.00 -0.00 ... 0.00 0.00 0.00]
[0.00 -0.00 1.00 ... 0.00 -0.00 0.00]
...
[-0.00 0.00 0.00 ... 1.00 0.00 0.00]
[0.00 0.00 -0.00 ... 0.00 1.00 0.00]
[-0.00 0.00 0.00 ... 0.00 0.00 1.00]]

This is A*1:  [[-2.00]
[-0.00]
[0.00]
[0.00]
[1.00]
[1.00]
[0.00]
[-1.00]
[0.00]
[-1.00]
[1.00]
[-1.00]
[-1.00]
[1.00]
[0.00]
.
.
.

What am I missing?