Construction Moore–Penrose inverse
1 construction
1.1 rank decomposition
1.2 qr method
1.3 singular value decomposition (svd)
1.4 block matrices
1.5 iterative method of ben-israel , cohen
1.6 updating pseudoinverse
1.7 software libraries
construction
rank decomposition
let
r
≤
min
(
m
,
n
)
{\displaystyle r\leq \min(m,n)}
denote rank of
a
∈
m
(
m
,
n
;
k
)
{\displaystyle a\in \mathrm {m} (m,n;k)\,\!}
.
a
{\displaystyle a\,\!}
can (rank) decomposed
a
=
b
c
{\displaystyle a=bc\,\!}
b
∈
m
(
m
,
r
;
k
)
{\displaystyle b\in \mathrm {m} (m,r;k)\,\!}
,
c
∈
m
(
r
,
n
;
k
)
{\displaystyle c\in \mathrm {m} (r,n;k)\,\!}
of rank
r
{\displaystyle r}
.
a
+
=
c
+
b
+
=
c
∗
(
c
c
∗
)
−
1
(
b
∗
b
)
−
1
b
∗
{\displaystyle a^{+}=c^{+}b^{+}=c^{*}(cc^{*})^{-1}(b^{*}b)^{-1}b^{*}\,\!}
.
the qr method
for
k
=
r
{\displaystyle k=\mathbb {r} \,\!}
or
k
=
c
{\displaystyle k=\mathbb {c} \,\!}
computing product
a
a
∗
{\displaystyle aa^{*}}
or
a
∗
a
{\displaystyle a^{*}a}
, inverses explicitly source of numerical rounding errors , computational cost in practice. alternative approach using qr decomposition of
a
{\displaystyle a\,\!}
may used instead.
consider case when
a
{\displaystyle a\,\!}
of full column rank,
a
+
=
(
a
∗
a
)
−
1
a
∗
{\displaystyle a^{+}=(a^{*}a)^{-1}a^{*}\,\!}
. cholesky decomposition
a
∗
a
=
r
∗
r
{\displaystyle a^{*}a=r^{*}r\,\!}
,
r
{\displaystyle r\,\!}
upper triangular matrix, may used. multiplication inverse done solving system multiple right-hand sides,
a
+
=
(
a
∗
a
)
−
1
a
∗
⇔
(
a
∗
a
)
a
+
=
a
∗
⇔
r
∗
r
a
+
=
a
∗
{\displaystyle a^{+}=(a^{*}a)^{-1}a^{*}\quad \leftrightarrow \quad (a^{*}a)a^{+}=a^{*}\quad \leftrightarrow \quad r^{*}ra^{+}=a^{*}}
which may solved forward substitution followed substitution.
the cholesky decomposition may computed without forming
a
∗
a
{\displaystyle a^{*}a\,\!}
explicitly, alternatively using qr decomposition of
a
=
q
r
{\displaystyle a=qr\,\!}
,
q
{\displaystyle q\,\,\!}
has orthonormal columns,
q
∗
q
=
i
{\displaystyle q^{*}q=i}
, ,
r
{\displaystyle r\,\!}
upper triangular. then
a
∗
a
=
(
q
r
)
∗
(
q
r
)
=
r
∗
q
∗
q
r
=
r
∗
r
{\displaystyle a^{*}a\,=\,(qr)^{*}(qr)\,=\,r^{*}q^{*}qr\,=\,r^{*}r}
,
so r cholesky factor of
a
∗
a
{\displaystyle a^{*}a}
.
the case of full row rank treated using formula
a
+
=
a
∗
(
a
a
∗
)
−
1
{\displaystyle a^{+}=a^{*}(aa^{*})^{-1}\,\!}
, using similar argument, swapping roles of
a
{\displaystyle a}
,
a
∗
{\displaystyle a^{*}}
.
singular value decomposition (svd)
a computationally simple , accurate way compute pseudoinverse using singular value decomposition. if
a
=
u
Σ
v
∗
{\displaystyle a=u\sigma v^{*}}
singular value decomposition of a,
a
+
=
v
Σ
+
u
∗
{\displaystyle a^{+}=v\sigma ^{+}u^{*}}
. rectangular diagonal matrix such
Σ
{\displaystyle \sigma }
, pseudoinverse taking reciprocal of each non-zero element on diagonal, leaving zeros in place, , transposing matrix. in numerical computation, elements larger small tolerance taken nonzero, , others replaced zeros. example, in matlab, gnu octave, or numpy function pinv, tolerance taken t = ε⋅max(m,n)⋅max(Σ), ε machine epsilon.
the computational cost of method dominated cost of computing svd, several times higher matrix–matrix multiplication, if state-of-the art implementation (such of lapack) used.
the above procedure shows why taking pseudoinverse not continuous operation: if original matrix has singular value 0 (a diagonal entry of matrix
Σ
{\displaystyle \sigma }
above), modifying may turn 0 tiny positive number, thereby affecting pseudoinverse dramatically have take reciprocal of tiny number.
block matrices
optimized approaches exist calculating pseudoinverse of block structured matrices.
the iterative method of ben-israel , cohen
another method computing pseudoinverse (cf. drazin inverse) uses recursion
a
i
+
1
=
2
a
i
−
a
i
a
a
i
,
{\displaystyle a_{i+1}=2a_{i}-a_{i}aa_{i},\,}
which referred hyper-power sequence. recursion produces sequence converging quadratically pseudoinverse of
a
{\displaystyle a}
if started appropriate
a
0
{\displaystyle a_{0}}
satisfying
a
0
a
=
(
a
0
a
)
∗
{\displaystyle a_{0}a=(a_{0}a)^{*}}
. choice
a
0
=
α
a
∗
{\displaystyle a_{0}=\alpha a^{*}}
(where
0
<
α
<
2
/
σ
1
2
(
a
)
{\displaystyle 0<\alpha <2/\sigma _{1}^{2}(a)}
,
σ
1
(
a
)
{\displaystyle \sigma _{1}(a)}
denoting largest singular value of
a
{\displaystyle a}
) has been argued not competitive method using svd mentioned above, because moderately ill-conditioned matrices takes long time before
a
i
{\displaystyle a_{i}}
enters region of quadratic convergence. however, if started
a
0
{\displaystyle a_{0}}
close moore–penrose inverse ,
a
0
a
=
(
a
0
a
)
∗
{\displaystyle a_{0}a=(a_{0}a)^{*}}
, example
a
0
:=
(
a
∗
a
+
δ
i
)
−
1
a
∗
{\displaystyle a_{0}:=(a^{*}a+\delta i)^{-1}a^{*}}
, convergence fast (quadratic).
updating pseudoinverse
for cases has full row or column rank, , inverse of correlation matrix (
a
a
∗
{\displaystyle aa^{*}}
full row rank or
a
∗
a
{\displaystyle a^{*}a}
full column rank) known, pseudoinverse matrices related
a
{\displaystyle a}
can computed applying sherman–morrison–woodbury formula update inverse of correlation matrix, may need less work. in particular, if related matrix differs original 1 changed, added or deleted row or column, incremental algorithms exist exploit relationship.
similarly, possible update cholesky factor when row or column added, without creating inverse of correlation matrix explicitly. however, updating pseudoinverse in general rank-deficient case more complicated.
software libraries
the python package numpy provides pseudoinverse calculation through functions matrix.i , linalg.pinv; pinv uses svd-based algorithm. scipy adds function scipy.linalg.pinv uses least-squares solver. high quality implementations of svd, qr, , substitution available in standard libraries, such lapack. writing 1 s own implementation of svd major programming project requires significant numerical expertise. in special circumstances, such parallel computing or embedded computing, however, alternative implementations qr or use of explicit inverse might preferable, , custom implementations may unavoidable.
the mass package r provides calculation of moore–penrose inverse through ginv function. ginv function calculates pseudoinverse using singular value decomposition provided svd function in base r package.
the octave programming language provides pseudoinverse through standard package function pinv pseudo_inverse() method.
Comments
Post a Comment