求解满足Ax=b的无唯一解酉矩阵A的Python实现问题

我需要求解形如$Ax=b$的线性方程组，其中A是未知的酉方阵，向量x和b已知。在物理问题中，A用于将随机单位向量x旋转为形式为$[1,0,0...0]$的b。该问题无唯一解，只需要其中一个可行解，且方案要能扩展至高维。

之前尝试过构造普通矩阵的方法，但得到的A不是酉矩阵；用cvxpy求解时，约束部分报错：LinAlgError：给定的是0维数组，数组必须至少是二维的，代码和报错信息如下：

测试代码 #

import numpy as np
import cvxpy as cp
n = 3
alpha=0.01
np.random.seed(1)
b= [1,0,0]
x = np.random.randn(n)+1j*np.random.randn(n)
x=x/np.linalg.norm(x)


#Construct the problem.
A = cp.Variable((n,n),complex=True)
objective = cp.Minimize(cp.sum_squares(A @ x - b))
constraints = [A.H==np.linalg.inv(A)]
prob = cp.Problem(objective, constraints)

# The optimal objective value is returned by `prob.solve()`.
result = prob.solve()
print(A.value)

报错信息 #

Traceback (most recent call last):
    File "/Users/sreerampg/untitled1.py", line 22, in <module>
        constraints = [A.H==np.linalg.inv(A)]
    File "<__array_function__ internals>", line 6, in inv
    File "/Applications/anaconda3/lib/python3.7/site-packages/numpy/linalg/linalg.py", line 539, in inv
        _assert_stacked_2d(a)
    File "/Applications/anaconda3/lib/python3.7/site-packages/numpy/linalg/linalg.py", line 197, in _assert_stacked_2d     
        'at least two-dimensional' % a.ndim)
LinAlgError：给定的是0维数组，数组必须至少是二维的

1. cvxpy约束错误原因 #

报错是因为np.linalg.inv(A)里的A是cvxpy的Variable对象，不是numpy数组，numpy无法对它求逆。酉矩阵的正确约束应该是A的共轭转置等于其逆，在cvxpy中要写成矩阵等式约束：A.H @ A == cp.eye(n)（等价于$A^H A = I$），或者直接用cvxpy内置的cp.is_unitary(A)约束（需cvxpy版本支持）。

2. 修复后的cvxpy代码 #

import numpy as np
import cvxpy as cp

n = 3
np.random.seed(1)
# 转成复数数组避免类型问题
b = np.array([1, 0, 0], dtype=np.complex128)
x = np.random.randn(n) + 1j * np.random.randn(n)
x = x / np.linalg.norm(x)

# 构造问题
A = cp.Variable((n, n), complex=True)
objective = cp.Minimize(cp.sum_squares(A @ x - b))
# 酉矩阵约束：共轭转置乘自身等于单位矩阵
constraints = [A.H @ A == cp.eye(n)]
prob = cp.Problem(objective, constraints)

# 使用支持复矩阵的求解器（如SCS）
result = prob.solve(solver=cp.SCS)
print("最优解矩阵A：")
print(A.value)

# 验证结果
print("\n验证Ax：", np.round(A.value @ x, 6))
print("验证A是否酉：", np.allclose(A.value.conj().T @ A.value, np.eye(n)))

3. 更高效的构造方法（无需优化） #

因为需求是把单位向量x映射到标准基向量$e_1=[1,0,...0]$，可以直接用Householder变换构造酉矩阵，这比凸优化更高效，且能直接扩展到任意高维：

原理 #

Householder变换可构造酉矩阵H，使得$Hx = e_1$，步骤如下：

1. 计算向量$v = x - e_1$
2. 若v的范数接近0（说明x已是$e_1$），则H为单位矩阵
3. 否则，构造$H = I - 2 \frac{vv^H}{|v|2}$，其中$v^H$是v的共轭转置

代码实现 #

import numpy as np

def householder_unitary(x):
    n = len(x)
    e1 = np.zeros(n, dtype=np.complex128)
    e1[0] = 1.0
    v = x - e1
    norm_v = np.linalg.norm(v)
    # 避免数值不稳定
    if norm_v < 1e-10:
        return np.eye(n, dtype=np.complex128)
    v = v / norm_v
    # 构造Householder酉矩阵
    H = np.eye(n, dtype=np.complex128) - 2 * np.outer(v, v.conj().T)
    return H

# 测试
n = 3
np.random.seed(1)
x = np.random.randn(n) + 1j * np.random.randn(n)
x = x / np.linalg.norm(x)

A = householder_unitary(x)
print("构造的酉矩阵A：")
print(A)

# 验证结果
print("\nAx结果：", np.round(A @ x, 6))
print("是否酉矩阵：", np.allclose(A.conj().T @ A, np.eye(n)))

该方法优势：

• 计算速度快，O(n²)复杂度，适配高维场景
• 直接构造，无需迭代优化
• 严格满足酉矩阵条件，数值误差范围内精确满足$Ax=e_1$

内容的提问来源于stack exchange，提问作者sreeram pg