Python和Matlab之间的正交、均值不变矩阵生成

Question

我需要一个程序来生成一个随机正交矩阵，它还保留矩阵列的平均值和协方差。换句话说，矩阵应具有任意正交矩阵A*A.T=I的性质，并且还具有A*1n=1n的均值保持性质，其中1n是大小为1xn的标量1的列。实现这一点的技术非常简单：

1. 生成一个随机矩阵
2. 连接一列1作为第一列
3. 对矩阵执行Gram-Schmidt正交归一化
4. 对另一个矩阵执行步骤1、2和3
5. 将这两个矩阵相乘

我在Matlab中有一个程序成功地实现了这一点：

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

上面的代码确实满足了我需要的两个属性

我想用Python编写类似的代码。我尝试了以下方法：

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(A@A.T))
print('This is A*1: ', np.around(A@onesVector))

尽管此Python代码满足第1个属性并给出了标识矩阵，但它不满足第2个属性，如下所示：

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]
.
.
.

我错过了什么