3.9 极小范数最小二乘广义逆

醒过来摸鱼

已于 2023-05-19 14:26:11 修改

阅读量1.6k

点赞数

分类专栏：线性代数【更新中】文章标签：线性代数矩阵算法

于 2023-02-07 21:50:51 首次发布

本文链接：https://blog.csdn.net/m0_66201040/article/details/128922340

版权

线性代数【更新中】专栏收录该内容

71 篇文章

订阅专栏

文章介绍了矩阵的广义逆，又称伪逆或穆尔-彭罗斯逆，其符号和性质，并提供了三种计算方法：广义逆法、满秩分解法和Greville迭代法。广义逆在解决最小二乘问题中有重要应用，且存在唯一性。文中还给出了相应的Python代码示例。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

定义

现在轮到最后一个广义逆了，这个广义逆名字特别多，最长的名字是极小范数最小二乘广义逆least square, minimum norm generalized inverse,短一点的名字是伪逆pseudoinverse，也叫穆尔-彭罗斯逆Moore-Penrose Inverse。它一般叫加号逆，符号为 $A^+$ ，也叫一二三四逆，符号为 $A^{\{1,2,3,4\}}$ ，因为符合所有的四个穆尔-彭罗斯条件。但是要注意这个加号和共轭转置很像，在西方的教材里，共轭转置的符号有时候是一个大宝剑 $A^{\dagger}$ ，其实是匕首符号，注意长短区分就好.
它的定义，其实我在别的文章里写过一遍了，就是 $G$ 如果是 $A$ 的加号逆，那么以下四个条件全部满足：
$AGA=A\\ GAG=G\\ (AG)^T=AG\\ (GA)^T=GA$
用加号逆求的 $A\bold x=\bold b$ 的解 $A^+\bold b$ ，是具有最小范数和最小二乘的解。关键是这样的逆存在吗？实际上，它是唯一存在的，过程我就不证明了。它的求解方法很多。我简要说三种算法。

广义逆法

广义逆法就是使用以下公式计算：
$A^+=A_m^-AA_l^-\\ =A^T(AA^T)^-A(A^TA)^-A^T$
其代码是五个矩阵相乘：

    # 极小范数最小二乘广义逆
    def min_norm_gls_inverse(self):
        a_t = self.transpose_matrix()
        left = (self * a_t).one_inverse()
        right = (a_t * self).one_inverse()
        return a_t * left * self * right * a_t

这样求得的逆，我举个例子：
$A=\begin{pmatrix}2 & 1 & 2\\ 2 & -2 & 2\\ 1 & 3 & 3\\ \end{pmatrix}\\ A^+=\begin{pmatrix}1 & -0.25 & -0.5\\ 0.333 & -0.333 & 0\\ -0.667 & 0.417 & 0.5\\ \end{pmatrix}$

满秩分解法

五个矩阵相乘，计算两个广义逆，想想还是挺可怕的。接下来介绍的算法，只需要四个矩阵相乘，但是还是要求两个逆矩阵。这种算法是利用矩阵的满秩分解来做，公式如下：
$A=BC\\ A^+=C^T(CC^T)^{-1}(B^TB)^{-1}B^T$
代码也不长：

    # 满秩分解法
    def moore_penrose_inverse(self):
        b, c = self.full_rank_factorization()
        ct = c.transpose_matrix()
        cct = c * ct
        bt = b.transpose_matrix()
        btb = bt * b
        return ct * cct.inverse() * btb.inverse() * bt

Greville迭代法

这个算法最早是卸载Greville于1960年发表的论文Some Applications of the Pseudoinverse of a Matrix里。这是一种按列迭代的算法，先计算第一列的加号逆，再计算前 $k$ 列组成的矩阵的加号逆，最后得到矩阵的加号逆。所以Greville算法的核心就是子矩阵递推关系。假设由前 $k$ 列组成的子矩阵为 $A_{k}$ ，其递推算法如下：
$d_k=A^+_{k-1}a_k\\ c_k=a_k-A_{k-1}d_k\\ b_k=\begin{cases} c_k^+, c_k\ne0\\ \bold (1+d_k^Td_k)^{-1}d_k^TA^+_{K-1}, c_k=0 \end{cases}\\ A_k^+=\begin{pmatrix}A^+_{k-1}-d_kb_k\\b_k\end{pmatrix}$
但是阅读Greville于1960年发表的论文时要注意，那个时代，匕首符号代表加号逆，不要认为是共轭转置。好了，代码我现在就贴出来：


    # greville算法
    def greville(self):
        vector_0 = self.__vectors[0]
        vector_0_inverse = Matrix([vector_0]).min_norm_gls_inverse()
        a_k_inverses = [vector_0_inverse]
        m = len(self.__vectors[0])
        n = len(self.__vectors)
        for k in range(1, n):
            d_k = a_k_inverses[k-1] * Matrix([self.__vectors[k]])
            a_k_1_matrix = self.sub_matrix(0, m, 0, k)
            c_k = matrix_utils.sub(self.__vectors[k], (a_k_1_matrix * d_k).__vectors[0])
            bk = None
            if matrix_utils.is_zero_vector(c_k):
                d_k_t = d_k.transpose_matrix()
                coefficient = 1 + (d_k_t * d_k).__vectors[0][0]
                bk = d_k_t* a_k_inverses[k-1]* coefficient
            else:
                bk = Matrix([c_k]).min_norm_gls_inverse()
            # 两个矩阵拼接起来
            up = a_k_inverses[k - 1] - (d_k * bk)
            a_k_inverse_array = [column + [bk.__vectors[i][0]] for i, column in enumerate(up.__vectors)]
            a_k_inverses.append(Matrix(a_k_inverse_array))
        return a_k_inverses[n-1]