最优化方法Python计算：构造正定矩阵的共轭向量

定义1 设$\mathbf{G}\in {\text{R}}^{n\times n}$为一正定矩阵。非零向量组${\mathbf{d}}_1,{\mathbf{d}}_2,\cdots ,{\mathbf{d}}_m\in {\text{R}}^n$（$m\le n$）满足
$${\mathbf{d}}_i^{\top }\mathbf{G}{\mathbf{d}}_j=0,1\le i\ne j\le m$$
称${\mathbf{d}}_1,{\mathbf{d}}_2,\cdots ,{\mathbf{d}}_m$关于$\mathbf{G}$共轭。
定理1 设无关向量组${\mathbf{\alpha }}_1,{\mathbf{\alpha }}_2,\cdots ,{\mathbf{\alpha }}_m\in {\text{R}}^n$，$\mathbf{G}\in {\text{R}}^{n\times n}$为一正定阵。令，
$${\mathbf{d}}_k=\{\begin{array}{ll}{\mathbf{\alpha }}_1& k=1\\ {\mathbf{\alpha }}_k-\sum _{i=1}^{k-1}\frac{{\mathbf{\alpha }}_k^{\top }\mathbf{G}{\mathbf{d}}_i}{{\mathbf{d}}_i^{\top }\mathbf{G}{\mathbf{d}}_i}{\mathbf{d}}_i& 1<k\le m\end{array}$$
则${\mathbf{d}}_1,{\mathbf{d}}_2,\cdots ,{\mathbf{d}}_m$关于$\mathbf{G}$共轭。
利用定理1，下列Python函数将${\text{R}}^n$中的自然基 { e 1 , e 2 , ⋯   , e n } \{\boldsymbol{e}_1,\boldsymbol{e}_2,\cdots,\boldsymbol{e}_n\} {e1​,e2​,⋯,en​}（线性无关）转换为关于指定正定矩阵$\mathbf{G}$共轭的向量组。

import numpy as np                                              #导入numpy
def conjugation(G):                                             #G为正定阵
    n,_=G.shape                                                 #G的行数
    A=np.eye(n)                                                 #初始线性无关组
    for k in range(1,n):                                        #依次构造G的共轭向量d2,...
        for i in range(k):                                      #按公式(3.3)构造dk
            numerator=np.matmul(np.matmul(A[:,k].T,G),A[:,i])   #第i项分子
            denominator=np.matmul(np.matmul(A[:,i].T,G),A[:,i]) #第i项分母
            A[:,k]-=numerator/denominator*A[:,i]                #减去第i项
    return A

程序的第2~10行定义函数conjugation，参数G表示正定阵$\mathbf{G}$，计算关于$\mathbf{G}$的共轭向量 { d 1 , d 2 , ⋯   , d n } \{\boldsymbol{d}_1,\boldsymbol{d}_2,\cdots,\boldsymbol{d}_n\} {d1​,d2​,⋯,dn​}。第3行读取G的行数n。第4行调用numpy的eye函数生成n阶单位阵${\mathbf{I}}_n$（每一列为一单位向量，所有列构成一线性无关组）赋予A。第5 ~9行的for循环，按定理1构造$\mathbf{G}$的共轭向量 { d k , k = 2 , ⋯   , n } \{\boldsymbol{d}_k,k=2,\cdots,n\} {dk​,k=2,⋯,n}（${\mathbf{d}}_1$即为A的第1列A[:,0]）。对表示第$k$个向量${\mathbf{d}}_k$的A的第k列A[:,k]，第6 ~9行的for循环按式(3.3)计算：第7行
$$\text{np.matmul(np.matmul(A[:,k].T,G),A[:,i])}$$
计算${\mathbf{\alpha }}_k^{\top }\mathbf{G}{\mathbf{d}}_i$赋予numerator，第8行
$$\text{np.matmul(np.matmul(A[:,i].T,G),A[:,i])}$$
计算${\mathbf{d}}_i^{\top }\mathbf{G}{\mathbf{d}}_i$赋予denominator，第9行计算
$$\text{numerator/denominator*A[:,i]}$$
计算$\frac{{\mathbf{\alpha }}_k^{\top }\mathbf{G}{\mathbf{d}}_i}{{\mathbf{d}}_i^{\top }\mathbf{G}{\mathbf{d}}_i}{\mathbf{d}}_i$并用其与A[:,k]（表示${\mathbf{\alpha }}_k$）相减。循环完成，按定理1，${\mathbf{d}}_k$的计算完成。执行函数完毕，返回存储在A中${\mathbf{G}}_{n\times n}$的共轭向量组 { d 1 , d 2 , ⋯   , d n } \{\boldsymbol{d}_1,\boldsymbol{d}_2,\cdots,\boldsymbol{d}_n\} {d1​,d2​,⋯,dn​}。conjugation函数可以用来构造正定矩阵的满秩共轭向量组。
例1 对正定阵$\mathbf{G}=\left(\begin{array}{ccc}3& 0& 1\\ 0& 4& 2\\ 1& 2& 3\end{array}\right)\in {\text{R}}^{3\times 3}$构造一组共轭相量 { d 1 , d 2 , d 3 } \{\boldsymbol{d}_1,\boldsymbol{d}_2,\boldsymbol{d}_3\} {d1​,d2​,d3​}。
解：下列代码完成计算。

import numpy as np                      #导入numpy
G=np.array([[3,0,1],                    #设置正定矩阵G
            [0,4,2],
            [1,2,3]],dtype='float')
print(conjugation(G))                   #计算并输出共轭向量组

程序很简单，借助注释信息不难理解。运行程序，输出

[[ 1.          0.         -0.33333333]
 [ 0.          1.         -0.5       ]
 [ 0.          0.          1.        ]]

即得到了关于$\mathbf{G}=\left(\begin{array}{ccc}3& 0& 1\\ 0& 4& 2\\ 1& 2& 3\end{array}\right)$的共轭向量组为${\mathbf{d}}_1=\left(\begin{array}{c}1\\ 0\\ 0\end{array}\right)$，${\mathbf{d}}_2=\left(\begin{array}{c}0\\ 1\\ 0\end{array}\right)$，${\mathbf{d}}_3=\left(\begin{array}{c}-\frac{1}{3}\\ -\frac{1}{2}\\ 1\end{array}\right)$。
写博不易，敬请支持：
如果阅读本文于您有所获，敬请点赞、评论、收藏，谢谢大家的支持！