矩阵最小二乘法推导

最小二乘法（Least Squares）用于求解线性方程组 的近似解，尤其是在方程组没有精确解时。其目标是最小化误差的平方和。只适用于非齐次方程解，对于齐次方程总有非平凡解

1. 问题描述

给定矩阵 和向量 ，我们希望找到向量 ，使得 ，即最小化以下目标函数：

其中，表示 向量的 2 范数，即欧几里得距离。

2. 目标函数展开

目标函数 是平方误差，可以展开如下：

展开后得到：

其中， 是常数项，和 无关。

3. 对 求导

为了找到使得 最小的 ，我们对 关于 求导数，并令其为 0。

求导得到：

令其等于 0：

4. 结果

通过上述推导，我们得到常见的 正规方程：

这是求解最小二乘问题的标准方程。

5. 解的存在与唯一性

• 存在性：正规方程有解的条件是 是可逆的。
• 唯一性：如果 可逆，则解是唯一的。通常，如果 的列满秩（即），则 是可逆的。

6. 解的形式

当 可逆时，最小二乘解可以表示为：

7. 总结

1. 最小二乘法通过最小化误差平方和来求解线性方程组的近似解。
2. 通过正规方程 可以求得最小二乘解。

SVD 求解最小二乘法

使用 奇异值分解（SVD）求解 最小二乘法 是一种常见的数值方法。它通过分解矩阵 ( A ) 来获得解，尤其在矩阵 ( A ) 不满秩或近似满秩时，使用 SVD 可以更稳定地求解。

1. 问题定义

给定矩阵 和向量 ，我们要解最小二乘问题：

这个问题的目标是找到一个向量 ，使得 尽量接近 。

2. 使用奇异值分解（SVD）

奇异值分解（SVD）将矩阵 ( A ) 分解为三个矩阵的乘积：

• 是一个 的正交矩阵（包含左奇异向量）；
• 是一个 的对角矩阵，其中对角线上的元素是奇异值；
• 是一个 的正交矩阵（包含右奇异向量）。

对于最小二乘问题，我们通过 SVD 进行分解后，可以得到解的表达式。

3. 最小二乘解的推导

3.1 使用 SVD 分解

我们首先用 SVD 将 ( A ) 分解：

要求解最小二乘问题 ，即：

由于 是正交矩阵，，我们可以将其作用到方程的两边：
得到：
这里 是新的向量，记作 ，于是方程变为：

3.2 解决方程

现在我们需要求解：

注意到 是对角矩阵，它包含了 的奇异值，记为 （其中 是 的秩）。我们可以将其表示为：

因此，我们有：

对每个 ( i )，有：

因此，得到：

所以最小二乘解 可以通过：

其中， 是 的伪逆，定义为：

即对于每个非零奇异值 ，其倒数成为 中相应位置的元素。

4. 最小二乘解的表达式

最终，最小二乘解 为:

其中，, , 是矩阵 的 SVD 分解结果， 是 的伪逆。

5. 数值解法

通过 SVD 来求解最小二乘法可以有效地处理奇异矩阵或者近似奇异矩阵的问题。SVD 是一种非常稳定的数值方法，尤其适用于解决欠定或超定线性方程组。

在实际应用中，可以使用 Python 的 numpy 或 scipy 库中的 svd 函数来进行 SVD 分解，然后计算最小二乘解。

示例代码（Python）：

import numpy as np

# 示例矩阵 A 和向量 b
A = np.array([[1, 2], [3, 4], [5, 6]])
b = np.array([7, 8, 9])

# 计算 SVD
U, Sigma, Vt = np.linalg.svd(A)

# 构造 Sigma 的伪逆
Sigma_plus = np.zeros((A.shape[1], A.shape[0]))  # 创建与 A.T 尺寸相同的矩阵
np.fill_diagonal(Sigma_plus, 1 / Sigma)          # 填充对角线

# 计算最小二乘解 x
x = Vt.T @ Sigma_plus @ U.T @ b

print("最小二乘解 x:")
print(x)

齐次矩阵最小二乘法

齐次矩阵最小二乘法是用来求解形如 的齐次方程组的问题。这里 是一个已知的矩阵，而 是待求解的向量。

1. 问题定义

给定矩阵 ，我们需要找到非零向量 ，使得：

直接解此方程可能存在平凡解（即 ）。为了避免平凡解，我们将问题重新表述为约束优化问题：

这里添加了 的约束，用于排除平凡解。

2. 使用奇异值分解（SVD）

2.1 矩阵 的 SVD 分解

对矩阵 进行奇异值分解：

其中：

• 是正交矩阵，包含左奇异向量；
• 是对角矩阵，包含奇异值 ；
• 是正交矩阵，包含右奇异向量。

2.2 最优解的确定

在约束优化问题中，目标是最小化 。将 替换为其 SVD 分解后：

由于 是正交矩阵，不改变向量的范数，因此有：

令 ，因为 也是正交矩阵，有 。约束条件 等价于 。

于是问题变为：

矩阵 的形式为：

其中 ，其余部分为零。

为了最小化 ， 应对应于 中最小奇异值的位置。因此，最优解为：

其中 是 的最后一列（对应于最小奇异值的右奇异向量）。

3. 结果总结

对于 的齐次问题，非平凡解可以通过 SVD 求解：

1. 计算 的奇异值分解 ；
2. 解 为 中对应最小奇异值的列向量。

4. 示例代码

以下是用 Python 实现的示例代码：

import numpy as np

# 示例矩阵 A
A = np.array([[1, 2, 3],
              [4, 5, 6],
              [7, 8, 9]])

# 计算 SVD
U, Sigma, Vt = np.linalg.svd(A)

# 最小奇异值对应的右奇异向量
h = Vt[-1, :]  # V 的最后一列是 Vt 的最后一行

print("非平凡解 h:")
print(h)

为最小奇异值对应右奇异向量的数学推导

给定矩阵 ，我们希望求解线性系统 的非零解。由线性代数理论，该问题可转化为寻找 使得：

1. 将问题转化为

设，代入奇异值分解 ，有：

由于 是正交矩阵，其列向量为正交基，满足 ，因此对齐次方程的求解不会影响解的性质。等式两边左乘 ，化简得：

令，有：

2. 的性质

根据 的结构， 等价于 的非零解只可能出现在最小奇异值对应的方向上，即最后一列的奇异向量方向。此时：

回代 可得：

结合 的取值， 是矩阵 的最后一列 ，即：

Homography 矩阵求解

Homography 矩阵（单应性矩阵）描述了两个平面之间的点到点的投影关系，是计算图像变换（如透视变换）的核心工具。设两幅图像中的对应点为 和 ，它们满足以下关系：

其中 是 的单应性矩阵，包含 8 个自由度（因比例因子可省略一个自由度）。

1. 问题定义

将齐次坐标展开后，得以下约束关系：

移项并整理为线性方程：

每对对应点 提供两条约束，求解 需要至少 4 对点（共 8 条约束）。

定义参数向量 ，上述方程可以写成矩阵形式：

其中 是由点坐标构成的 矩阵。

2. 使用 SVD 求解

要解方程 ，采用齐次最小二乘法：

2.1 奇异值分解

对矩阵 进行奇异值分解（SVD）：

其中：

• 是正交矩阵，包含左奇异向量；
• 是对角矩阵，包含奇异值；
• 是正交矩阵，包含右奇异向量。

2.2 最优解

取 的最后一列（对应最小奇异值的右奇异向量）作为 ，即：

将其重新排列为 矩阵，得到 。

3. 归一化处理

由于 是齐次矩阵，不受比例因子的影响，计算结果常需要归一化处理。将 的元素整体除以 ，使得：

4. 实现代码

以下是用 Python 实现 Homography 矩阵求解的代码：

import numpy as np

def compute_homography(points_src, points_dst):
    """
    计算单应性矩阵 H，使得 points_dst ≈ H * points_src。
    
    :param points_src: 源图像点坐标 (Nx2)
    :param points_dst: 目标图像点坐标 (Nx2)
    :return: 单应性矩阵 H (3x3)
    """
    assert points_src.shape == points_dst.shape, "点集尺寸不一致"
    assert points_src.shape[0] >= 4, "至少需要4对点"
    
    n = points_src.shape[0]
    A = []
    
    for i in range(n):
        x, y = points_src[i]
        x_prime, y_prime = points_dst[i]
        A.append([-x, -y, -1, 0, 0, 0, x*x_prime, y*x_prime, x_prime])
        A.append([0, 0, 0, -x, -y, -1, x*y_prime, y*y_prime, y_prime])
    
    A = np.array(A)
    
    # SVD 分解
    U, Sigma, Vt = np.linalg.svd(A)
    
    # 最小奇异值对应的右奇异向量
    H = Vt[-1, :].reshape((3, 3))
    
    # 归一化，使得 H[2, 2] = 1
    H /= H[-1, -1]
    return H

# 示例点
points_src = np.array([[0, 0], [1, 0], [1, 1], [0, 1]])
points_dst = np.array([[0, 0], [2, 0], [2, 2], [0, 2]])

H = compute_homography(points_src, points_dst)
print("单应性矩阵 H:")
print(H)

图片压缩

import cv2
import numpy as np

top_k = 360  # 保留的奇异值数量

# 读取彩色图像
image = cv2.imread("test.jpg")
print(image.shape)  # 输出图像的形状

# 分离颜色通道 (BGR)
blue_channel = image[:, :, 0]
green_channel = image[:, :, 1]
red_channel = image[:, :, 2]

# 对每个通道进行SVD
U_b, Sigma_b, Vt_b = np.linalg.svd(blue_channel, full_matrices=False)
U_g, Sigma_g, Vt_g = np.linalg.svd(green_channel, full_matrices=False)
U_r, Sigma_r, Vt_r = np.linalg.svd(red_channel, full_matrices=False)

# 保留前 top_k 个奇异值，并重构每个通道
Sigma_b = np.diag(Sigma_b[:top_k])
Sigma_g = np.diag(Sigma_g[:top_k])
Sigma_r = np.diag(Sigma_r[:top_k])

blue_reconstructed = np.dot(U_b[:, :top_k], np.dot(Sigma_b, Vt_b[:top_k, :]))
green_reconstructed = np.dot(U_g[:, :top_k], np.dot(Sigma_g, Vt_g[:top_k, :]))
red_reconstructed = np.dot(U_r[:, :top_k], np.dot(Sigma_r, Vt_r[:top_k, :]))

# 将重构的通道合并为一张彩色图像
new_image = np.stack([blue_reconstructed, green_reconstructed, red_reconstructed], axis=-1)

# 显示重构的彩色图像
cv2.imshow("Reconstructed Image", new_image.astype(np.uint8))
cv2.waitKey(0)
cv2.destroyAllWindows()