求解具有未知数约束的线性方程组

Question

我想求解未知向量 $h$ 的线性方程组 $y=Uh$，其中我对 $h$ 的某些元素有一些约束。

让我们考虑一个小例子：$y$ 是大小为 $N$ 的向量，$U$ 是大小为 $N imes 9$ 的矩阵，$h$ 是大小为 $9$ 的向量，其中 $N>9 $。在这里，我对向量 $h$ 有一些限制。如果前三个元素是 $a、b、c$，那么接下来的三个元素应该是 $pa^2、pb^2、pc^2$，其中 $p$ 是非零缩放因子。最后三个元素应为 $2pab、2pac、2pbc$。如何使用未知向量的另一个参数对参数使用此类约束？有人可以帮我解决这个问题吗？

Answer 1

这不是一个线性系统；这是一个非线性系统。此外，您几乎肯定必须采用较弱的 solve 定义，因为系统是超定的，因此只能在最小二乘意义上求解。只有三个未知数（a、b、c），且 N>9 方程。举例来说，

import numpy as np
from scipy.optimize import least_squares


def make_h(abc: np.ndarray, p: float) -> np.ndarray:
    a, b, c = abc

    # Here, I have few constraints on the vector $h$.
    return np.concatenate((
        # If the first three elements are $a, b, c$,
        abc,
        # then the next three elements should be $pa^2, pb^2, pc^2$.
        p*abc**2,
        # The last three elements should then be $2pab, 2pac, 2pbc$.
        (
            2*p*a*b,
            2*p*a*c,
            2*p*b*c,
        ),
    ))


def residuals(
    abc: np.ndarray, p: float, U: np.ndarray, y: np.ndarray,
) -> np.ndarray:
    h = make_h(abc, p)

    # I would like to solve a system of linear equations $y=Uh$ for an unknown vector $h$.
    yhat = U @ h
    return yhat - y


def solve(
    p: float, U: np.ndarray, y: np.ndarray,
) -> tuple[np.ndarray, np.ndarray]:
    result = least_squares(
        fun=residuals,
        args=(p, U, y),
        x0=(1, 1, 1),
    )
    if not result.success:
        raise ValueError(result.message)

    return make_h(result.x, p), result.x


def main() -> None:
    rand = np.random.default_rng(seed=0)
    N = 20                        # $N>9$.
    U = rand.random(size=(N, 9))  # $U$ is a matrix of size $N \times 9$
    p = rand.random()             # $p$ is a non-zero scaling factor.

    # Make some fake data that stand a chance of being solvable
    hidden_abc = rand.random(size=3)
    hidden_h = make_h(hidden_abc, p)
    noise = rand.normal(size=N, scale=0.05)

    # $y$ is a vector of size $N$
    y = noise + U@hidden_h

    h, abc = solve(p, U, y)
    print('h =', h)
    print('abc =', abc)


if __name__ == '__main__':
    main()