概要

　有限要素法（FEM）や粒子法を使い、物理量を求める際、連立一次方程式に帰着して解くことがあります。

　連立一次方程式に、拘束条件を入れてラグランジュ未定乗数法で解く場合、支配方程式から得られる剛性行列にラグランジュ未定乗数に関わる行列が追加されるので大規模行列になります。剛性行列の自由度を増やさず、代数的に解く方法を紹介します。

論文は以下を引用しています。

ラグランジュ未定乗数法

解きたい支配方程式は、

\begin{align}
\sum_jS_{ij}x_j=f_i
\end{align}

であり（$S$は正定値行列？）、拘束条件は

\begin{align}
\sum_jC_{ij}x_j=g_i
\end{align}

とする。目的関数を、拘束条件の下で最小化することを考える。ラグランジュ未定乗数を使い、目的関数$F(x)$を

\begin{align}
F(x) = \frac{1}{2}\sum_{i,j}S_{ij}x_ix_j - \sum_i f_ix_i + \sum_i \lambda_i \left(\sum_j C_{ij}x_j -g_i \right)
\end{align}

とする。目的関数$F(x)$が最小となる$x$と$\lambda$は、

\begin{align}
\frac{\partial F(x)}{\partial x_k} &= \frac{1}{2}\left(\sum_{j}S_{kj}x_j +  \sum_{i}S_{ik}x_i  \right) - f_k + \sum_i \lambda_i  C_{ik} =0 \\
\frac{\partial F(x)}{\partial \lambda_k} &= \sum_j C_{kj}x_j -g_k  =0
\end{align}

を満たす。$S$を対称行列とすれば、

\begin{pmatrix}
S & C^T \\
C & 0  \\
\end{pmatrix}
\begin{pmatrix}
x \\
\lambda \\
\end{pmatrix} = 
\begin{pmatrix}
f \\
g \\
\end{pmatrix}

を解けばよい。剛性行列$S$にラグランジュ未定乗数に関わる行列$C$が追加されるので大規模行列になる。この行列を直接解くことはせず、代数的に解く方法を以下説明する。

　定理1

　$S$を正方行列、$C$を列の大きさが$S$と同じ横長の行列とする。$P,Q,R$ を

\begin{align}
P&=C^T(CC^T)^{-1}C\\
Q&=I-P \\
R&=C^T(CC^T)^{-1}
\end{align}

とする。2つの条件
1．$C$がフルランク
2. $\mbox{Ker}(QSQ)\cap \mbox{Ker}(C) ={0}$
が成立するとき、$S_1$と$f_1$を

\begin{align}
S_1&=C^TC +QSQ　\\
f_1&=C^Tg+Q(f-SRg)
\end{align}

とおくと、
1．$S_1x=f_1$は一意に決まる。
2. $\lambda=R^T(f-Sx)$
が成立する。

　$P$は、射影行列であり、

\begin{align}
PP&=C^T(CC^T)^{-1}CC^T(CC^T)^{-1}C=P \\
QQ&=(I-P)(I-P)=I-P-P+P=Q
\end{align}

となる。

　$\mbox{Ker}(A)$は、零空間（核空間）と呼ばれ

\begin{align}
\mbox{Ker}(A)= \{x\in V ; Ax=0\}
\end{align}

$A$は線形作用素で、$A:V\rightarrow W$、$V,W$はベクトル空間である。また

\begin{align}
\mbox{Aが正則} \Leftrightarrow \mbox{ker}(A)=\{0\}
\end{align}

である。以下の補題2ではこの関係を使う。

補題2

$\mbox{Ker}(QSQ)\cap \mbox{Ker}(C) =\mbox{{0}}$ のとき、$CC^T+QSQ$は正則。

補題2の証明

　
$(C^TC+QSQ)x=0$とすると、左から$Q$を掛けて

\begin{align}
(QC^TC+Q^2SQ)x=0
\end{align}

となるが、

\begin{align}
QC^T = (I-C^T(CC^T)C)C^T =C^T-C^T=0
\end{align}

なので、

\begin{align}
&Q^2SQx=QSQx=0 \\
\therefore \ \ & x \in \mbox{Ker} (QSQ)
\end{align}

となる。

同様に$(C^TC+QSQ)x=0$とすると、左から$P$を掛けて

\begin{align}
(PC^TC+PQSQ)x=0
\end{align}

となるが、

\begin{align}
PQ=P-P=0
\end{align}

なので、$C$はフルランクであることを用いると

\begin{align}
&PC^TCx=C^T(CC^T)^{-1}CC^TCx=C^TCx = 0 \\
&Cx=0 \\
\therefore \ \ & x \in \mbox{Ker} (C)
\end{align}

となる。

　仮定より、$\mbox{Ker}(QSQ)\cap \mbox{Ker}(C) =\mbox{{0}}$より、$x=0$となる。よって$C^TC+QSQ$は正則。

定理1の証明

　以下の行列$A$を考える。

A=
\begin{pmatrix}
Q & C^T-QSR \\
R^T & 0  \\
\end{pmatrix}

この逆行列は、

A^{-1}=
\begin{pmatrix}
Q(I+SRR^T) & C^T \\
R^T & 0  \\
\end{pmatrix}

であることを確認する。まず

\begin{align}
QC^T&=0\\
PQ& =0 \\
RC&=C^T(CC^T)^{-1}C=P \\
C^TR^T&=C^T(CC^T)^{-1}C=P \\
R^TC^T&=(CR)^T=(CC^T(CC^T)^{-1})=1 \\
R^TQ&=(C^T(CC^T)^{-1})^T(I-P)\\
&=(CC^T)^{-T}C-(CC^T)^{-T}CC^T(CC^T)^{-1}C=0
\end{align}

であることを用いれば、

\begin{pmatrix}
Q(I+SRR^T) & C^T \\
R^T & 0  \\
\end{pmatrix}
\begin{pmatrix}
Q & C^T-QSR \\
R^T & 0  \\
\end{pmatrix}

の$(1,1)$の成分は、

\begin{align}
QQ+SRR^TQ+C^TR^T = Q+P=I
\end{align}

$(1,2)$の成分は、

\begin{align}
QC^T-Q^2SR+QSRR^TC^T-QSRR^TQSR=-QSR+QSR=0
\end{align}

$(2,1)$の成分は、

\begin{align}
R^TQ=0
\end{align}

$(2,2)$の成分は、

\begin{align}
R^TC^T-R^TQSR=I
\end{align}

となり、$A$は逆行列が存在する。

　次に、本来解きたい行列に左から$A$を掛けた

\begin{align}
\begin{pmatrix}
Q & C^T-QSR \\
R^T & 0  \\
\end{pmatrix}
\left\{
\begin{pmatrix}
S & C^T \\
C & 0  \\
\end{pmatrix}
\begin{pmatrix}
x \\
\lambda \\
\end{pmatrix} - 
\begin{pmatrix}
f \\
g \\
\end{pmatrix}
\right\}=0
\end{align}

を考える。$A$は逆行列が存在するので、上記の行列を解くことと、本来解きたい行列を解くことは同じである。この行列を計算していくと、

\begin{align}
\begin{pmatrix}
QS+C^TC-QSRC & QC^T \\
R^TC & R^TC^T  \\
\end{pmatrix}
\begin{pmatrix}
x \\
\lambda \\
\end{pmatrix} 
=
\begin{pmatrix}
Qf+C^Tg-QSRg \\
R^Tf \\
\end{pmatrix}
\end{align}

\begin{align}
\begin{pmatrix}
QS+C^TC-QS(I-Q) & 0 \\
R^TC & I  \\
\end{pmatrix}
\begin{pmatrix}
x \\
\lambda \\
\end{pmatrix} 
=
\begin{pmatrix}
C^Tg+Q(f-SRg) \\
R^Tf \\
\end{pmatrix}
\end{align}

\begin{align}
\begin{pmatrix}
C^TC+QSQ & 0 \\
R^TC & I  \\
\end{pmatrix}
\begin{pmatrix}
x \\
\lambda \\
\end{pmatrix} 
=
\begin{pmatrix}
C^Tg+Q(f-SRg) \\
R^Tf \\
\end{pmatrix}
\end{align}

となる。これは、

\begin{align}
S_1&=C^TC +QSQ　\\
f_1&=C^Tg+Q(f-SRg)
\end{align}

とおくと、$x,\lambda$は、

\begin{align}
&S_1x=f_1\\
&R^TSx+I\lambda=R^Tf \\
&\therefore \lambda = R^T(f-Sx)
\end{align}

となり、$\lambda = R^T(f-Sx)$が得られる。補題2より$S_1=C^TC+QSQ$は正則なので、$S_1x=f_1$は一意に解が求まる。

まとめ

　定理1より、

\begin{align}
&S_1x=f_1\\
&\lambda = R^T(f-Sx)
\end{align}

を解けばよいが、実際にほしい物理量は$x$なので、コードは以下のように書けばよいだろう（定理１で必要な２つの条件は考慮していないが）。
　
　定理1では、剛性行列$S$が、正定値対称行列であることは仮定されなくても証明される。

　$C$がフルランクではない場合、ムーア-ペンローズの擬似逆行列を使い計算できることが、論文に記載してある。機会があれば紹介したい。

S = 行列
f = ベクトル
C = Sと列の長さが同じ横長行列
g = ベクトル

CT=C.T

# CCT 逆行列
CCT = np.matmul(C,CT)
inv_CCT = inv(csc_matrix(CCT)).toarray()
# P = CT(CCT)^{-1}C
P = np.matmul(CT,np.matmul(inv_CCT,C))
# Q = 1 - P
Q =np.eye(P.shape[0]) -P
# R = CT(CCT)^{-1}
R = np.matmul(CT,inv_CCT)
# CT*g
CTg=np.matmul(CT,g)
# Rg
Rg = np.matmul(R,g)

S1 = np.matmul(CT,C) + np.matmul(Q,np.matmul(S,Q))
f1 = CTg + np.matmul(Q, f-np.matmul(S,Rg) )

Mat=lil_matrix(S1).tocsr()
print(Mat)
# 解
from scipy.sparse.linalg.dsolve import linsolve
x = linsolve.spsolve(Mat, f1)