関数singularvalue_decomposition マニュアル

(The documentation of function singularvalue_decomposition)
Last Update: 2021/12/6


◆機能・用途(Purpose)

正方行列または縦長の矩形行列の特異値分解を計算する。
Calculate the singular value decomposition of a square or vertically long matrix.

◆形式(Format)

#include <matrix/inverse.h>
inline struct matrix singularvalue_decomposition(const struct matrix G)


◆引数(Arguments)

G 特異値分解をしたい行列\(\myvector{G}\)。正方形または縦長のもののみ可。
The matrix, \(\myvector{G}\), for which the singular value decomposition is needed, which must be a square or vertically long matrix.


◆戻り値(Return value)

行列\(\myvector{G}\)の特異値分解を \[\begin{equation} \myvector{G}=\myvector{U} \myvector{\Lambda} \myvector{V}^T \label{eq.decomposition} \end{equation}\] としたとき、行列 \[\begin{equation} \begin{pmatrix} \myvector{\Lambda} & \myvector{U}^T \\ \myvector{V} & \myvector{0} \end{pmatrix} \label{eq.returnValue} \end{equation}\] を戻り値として返す。
A matrix given by eq. (\ref{eq.returnValue}) is returned, where the matrices \(\myvector{U}\), \(\myvector{\Lambda}\), and \(\myvector{V}\) are defined by the singular value decomposition (eq. \ref{eq.decomposition}).
ここで上付き添字の\(^T\)は行列の転置を表し、 \(\myvector{G}\)のサイズを\(N\times M\)としたとき \(\myvector{U}\), \(\myvector{\Lambda}\), \(\myvector{V}\)はそれぞれ \(N\times N\), \(N\times M\), \(M\times M\)の大きさの行列として 以下の式で定義される。 \[\begin{equation} \myvector{U}= \begin{pmatrix} \myvector{u_1} & \cdots & \myvector{u_N} \end{pmatrix} \label{eq.U} \end{equation}\] \[\begin{equation} \myvector{\Lambda}= \begin{pmatrix} \lambda_1 & 0 & \cdots & 0 \\ 0 & \ddots & \ddots & \vdots \\ \vdots & \ddots & \ddots & 0 \\ 0 & \cdots & 0 & \lambda_M \\ 0 & \cdots & \cdots & 0 \\ \vdots & & & \vdots \\ 0 & \cdots & \cdots & 0 \end{pmatrix} \label{eq.Lambda} \end{equation}\] \[\begin{equation} \myvector{V}= \begin{pmatrix} \myvector{v_1} & \cdots & \myvector{v_M} \end{pmatrix} \label{eq.V} \end{equation}\] ここで\(\lambda_i\), \(\myvector{u_i}\), \(\myvector{v_i}\) (\(i=1,\cdots,M\))は固有値問題 \[\begin{equation} \begin{pmatrix} \myvector{0} & \myvector{G} \\ \myvector{G}^T & \myvector{0} \end{pmatrix} \begin{pmatrix} \myvector{u_i} \\ \myvector{v_i} \end{pmatrix} =\lambda_i \begin{pmatrix} \myvector{u_i} \\ \myvector{v_i} \end{pmatrix} \label{eq.enlarged_G} \end{equation}\] に対する固有値, 固有ベクトルである。 (\ref{eq.enlarged_G})を成分で書き下すと \(\myvector{G}\myvector{v_i}=\lambda_i\myvector{u_i}\), \(\myvector{G}^T\myvector{u_i}=\lambda_i\myvector{u_i}\) となることから \[\begin{equation} \myvector{G}^T\myvector{G}\myvector{v_i} =\myvector{G}^T\lambda_i\myvector{u_i} =\lambda_i\myvector{G}^T\myvector{u_i} =\lambda_i^2\myvector{v_i} \label{eq.eigen.GTG} \end{equation}\] \[\begin{equation} \myvector{G}\myvector{G}^T\myvector{u_i} =\myvector{G}\lambda_i\myvector{v_i} =\lambda_i\myvector{G}\myvector{v_i} =\lambda_i^2\myvector{u_i} \label{eq.eigen.GGT} \end{equation}\] であるので、\(\myvector{v_i}\)は 行列\(\myvector{G}^T\myvector{G}\)の固有ベクトル、 \(\myvector{u_i}\)は行列\(\myvector{G}\myvector{G}^T\)の固有ベクトル、 \(\lambda_i^2\)はそれらの共通の固有値という見方もできる。 このような見方をすれば \(\myvector{u_1},\cdots,\myvector{u_M}\)が \(\myvector{G}\myvector{G}^T\)の\(N\)個の固有ベクトルのうちの 固有値\(\lambda_i^2\)に対応する\(M\)個であることから、 残りの\(N-M\)個の固有ベクトルとして \(\myvector{u_{M+1}},\cdots,\myvector{u_N}\)が自然に定義される。 なお\(\myvector{u_{M+1}},\cdots,\myvector{u_N}\)は 固有値0に対する固有ベクトルである。
Here, the superscript \(^T\) represents the transpose of a matrix, and \(\myvector{U}\), \({\Lambda}\), and \(\myvector{V}\) are \(N\times N\), \(N\times M\), and \(M\times M\) matrices defined by eqs (\ref{eq.U})-(\ref{eq.V}), respectively, where \(N\) and \(M\) are the numbers of rows and columns of \(\myvector{G}\). In this definition, \(\lambda_i\), \(\myvector{u_i}\), and \(\myvector{v_i}\) (\(i=1,\cdots,M\)) are eigenvalues and eigenvectors of an enlarged matrix given by eq. (\ref{eq.enlarged_G}). From this equation, we have \(\myvector{G}\myvector{v_i}=\lambda_i\myvector{u_i}\) and \(\myvector{G}^T\myvector{u_i}=\lambda_i\myvector{u_i}\), and thus (\ref{eq.eigen.GTG}) and (\ref{eq.eigen.GGT}), indicating that \(\myvector{v_i}\) and \(\myvector{u_i}\) for \(i=1,\cdots,M\) are eigenvectors of matrices \(\myvector{G}^T\myvector{G}\) and \(\myvector{G}\myvector{G}^T\), respectively, and \(\lambda_i^2\) is a common eigenvalue for them. On this viewpoint, \(\myvector{u_1},\cdots,\myvector{u_M}\) are the first \(M\) eigenvectors corresponding to eigenvalues \(\lambda_i^2\), so that \(\myvector{u_{M+1}},\cdots,\myvector{u_N}\) are defined naturally as the remaining \(N-M\) eigenvectors corresponding to eigenvalues of 0.

◆使用例(Example)

struct matrix G,U,UT,Lambda,V,VT,tmp;
tmp=singularvalue_decomposition(G);
Lambda=cut_matrix(tmp,1,G.rowmax,1,G.columnmax);
UT=cut_matrix(tmp,1,G.rowmax,G.columnmax+1,tmp.columnmax);
V=cut_matrix(tmp,G.rowmax+1,tmp.rowmax,1,G.columnmax);
U=matrix_transpose(UT);
VT=matrix_transpose(V);

この例ではU,Lambda,VTがそれぞれ \(\myvector{U}\),\(\myvector{\Lambda}\),\(\myvector{V}^T\)になる。
In this example, “U”, “Lambda”, and “VT” represent \(\myvector{U}\), \(\myvector{\Lambda}\), and \(\myvector{V}^T\), respectively.

◆計算方法(Computation method)

この関数では以下の手順で特異値分解の計算を行っている。
In this function, the singular value decomposition is investigated by the following steps:
  1. 計算を安定させるために\(\myvector{G}\)の全要素を絶対値最大値で規格化する。
    Normalize all the components of \(\myvector{G}\) by the absolute maximum to stabilize the calculation.
  2. Householder変換を用いて\(\myvector{G}\)を2重対角行列にする (関数double_diagonal_householder (matrix/householder.h)使用)。
    The matrix \(\myvector{G}\) is converted to a double diagonal matrix using Householder conversion by function double_diagonal_householder (matrix/householder.h).
  3. 得られた2重対角行列をGivens変換を用いて対角行列にする (関数double_diagonal_2_diagonal (matrix/givens.h)使用)。
    The double diagonal matrix is converted to a diagonal matrix using Givens conversion by function double_diagonal_2_diagonal (matrix/givens.h).
  4. 特異値を左上から右下に向かって降順になるようにソートする。
    Sort the singular values to a descending order from left top to right bottom.
  5. 特異値行列\(\myvector{\Lambda}\)を 規格化前の\(\myvector{G}\)に対するものに直す。
    Correct the singular value matrix \(\myvector{\Lambda}\) to what would be obtained if \(\myvector{G}\) was not normalized.
  6. 各特異値\(\lambda_i\) (\(i=1,\cdots,M\))の符号をチェックし、 もし負であれば\(\lambda_i\)と対応する固有ベクトル\(\myvector{u_i}\) の符号を\(-1\)倍する。
    Sign of each singular value \(\lambda_i\) (\(i=1,\cdots,M\)) is checked. If the sign is negative, \(\lambda_i\) and the corresponding eigenvector \(\myvector{u_i}\) are multiplied with \(-1\).

なお\(\lambda_i\)と\(\myvector{u_i}\)を同時に\(-1\)倍すれば (\ref{eq.decomposition})式の関係自体は保存される。 このことは(\ref{eq.decomposition})式を \[\begin{equation} \myvector{U}= \begin{pmatrix} \myvector{u_1} & \cdots & \myvector{u_N} \end{pmatrix} \begin{pmatrix} \lambda_1 & 0 & \cdots & 0 \\ 0 & \ddots & \ddots & \vdots \\ \vdots & \ddots & \ddots & 0 \\ 0 & \cdots & 0 & \lambda_M \\ 0 & \cdots & \cdots & 0 \\ \vdots & & & \vdots \\ 0 & \cdots & \cdots & 0 \end{pmatrix} \begin{pmatrix} \myvector{v_1}^T \\ \vdots \\ \myvector{v_M}^T \end{pmatrix} =\sum_{i=1}^M \lambda_i\myvector{u_i}\myvector{v_i^T} \label{eq.decomposition.component} \end{equation}\] と成分表記してみれば簡単に分かる。 このことから特異値, 特異ベクトルには\(-1\)倍の不確定性があり、 負の特異値が登場する場合があるのはこの不確定性が原因と思われる。 参考までに、固有値分解の場合には\(\myvector{u_i}=\myvector{v_i}\)なので \(\myvector{v_i}\)の符号を変えずに\(\lambda_i\)と\(\myvector{u_i}\)のみ \(-1\)倍するという操作は出来ず、よって\(-1\)倍の不確定は生じない。
Note that the relation of eq. (\ref{eq.decomposition}) is preserved by simultaneously multiplying \(-1\) to \(\lambda_i\) and \(\myvector{u_i}\). This is easily confirmed by looking at eq. (\ref{eq.decomposition.component}). Therefore the singular values and vectors have an ambiquity of the sign, which seems to be the cause for negative singular values. For comparison, the ambiguity does not occur in case of an eigenvalue decomposition, for which a relation \(\myvector{u_i}=\myvector{v_i}\) holds; awing to this relation, it is impossible to multiply \(-1\) to \(\lambda_i\) and \(\myvector{u_i}\) keeping the sign of \(\myvector{v_i}\).

◆検証(Validation)

この関数が正しい結果を返すことを確かめるには、 この関数で得られた\(\myvector{U}\), \(\myvector{\Lambda}\), \(\myvector{V}^T\)を用いて (\ref{eq.decomposition})式の \(\myvector{U}\myvector{\Lambda}\myvector{V}^T\)を計算し、 元の\(\myvector{G}\)と同じになることを確かめれば良い。 乱数を用いて作成した以下の行列を用いてこの検証を行った。
To check that this function returns a correct singular value decomposition, calculate \(\myvector{U}\myvector{\Lambda}\myvector{V}^T\) (eq. \ref{eq.decomposition}) using the matrices \(\myvector{U}\), \(\myvector{\Lambda}\), \(\myvector{V}^T\) obtained by this function, and compare the result with the given matrix \(\myvector{G}\). This examination was conducted using the matrix below which was generated by random values.
\[\begin{equation} \myvector{G}= \begin{pmatrix} -5.542348e-10 & 1.180734e-09 & -5.946389e-10 & -5.832139e-10 & -5.107910e-10 \\ -2.357166e-09 & 1.389110e-09 & 6.061486e-10 & -1.676399e-09 & -8.405894e-10 \\ 9.754172e-10 & -7.404720e-10 & 1.276538e-09 & -9.070130e-10 & -4.890220e-10 \\ 5.082556e-10 & 7.325038e-10 & 6.491890e-10 & 3.288509e-09 & 7.671914e-10 \\ -2.856718e-08 & 1.917195e-09 & -3.393252e-09 & -5.790534e-10 & 2.972070e-09 \\ 1.161411e-09 & 3.587793e-09 & -4.279629e-09 & -4.661626e-10 & -2.133547e-09 \\ -9.078414e-10 & 5.549452e-10 & 7.600898e-10 & 1.573012e-09 & -7.303892e-10 \end{pmatrix} \end{equation}\]
この\(\myvector{G}\)の特異値分解をこの関数で求めたところ、以下の結果となった。
The singular value decomposition of this \(\myvector{G}\) obtained by this function was as follows:
\[\begin{equation} \myvector{U}= \begin{pmatrix} -2.257041e-02 & -2.122645e-01 & -7.089979e-02 & 2.103495e-01 & -2.733569e-01 & -8.498114e-01 & 3.290246e-01 \\ -7.917271e-02 & -1.251547e-01 & -3.631409e-01 & 7.271012e-01 & -4.125188e-01 & 3.834229e-01 & 1.831903e-02 \\ 4.064394e-02 & 1.417409e-01 & -2.960944e-01 & 2.544087e-01 & 2.993850e-01 & -3.486574e-01 & -7.840080e-01 \\ 1.805492e-02 & 1.299719e-01 & 8.001442e-01 & 2.204498e-01 & -4.158117e-01 & -4.495841e-02 & -3.450098e-01 \\ -9.950663e-01 & 1.529178e-03 & 2.748692e-02 & -6.332450e-02 & 2.916873e-02 & -2.459690e-02 & -6.016154e-02 \\ 2.117112e-02 & -9.487043e-01 & 1.460286e-01 & -4.711964e-02 & 1.378046e-01 & 8.134367e-02 & -2.244111e-01 \\ -2.528293e-02 & 4.748310e-02 & 3.363322e-01 & 5.545491e-01 & 6.875554e-01 & 1.547192e-03 & 3.220670e-01 \end{pmatrix} \end{equation}\] \[\begin{equation} \myvector{\Lambda}= \begin{pmatrix} 2.913348e-08 & 0.000000e+00 & 0.000000e+00 & 0.000000e+00 & 0.000000e+00 \\ 0.000000e+00 & 6.370206e-09 & 0.000000e+00 & 0.000000e+00 & 0.000000e+00 \\ 0.000000e+00 & 0.000000e+00 & 4.191569e-09 & 0.000000e+00 & 0.000000e+00 \\ 0.000000e+00 & 0.000000e+00 & 0.000000e+00 & 2.367502e-09 & 0.000000e+00 \\ 0.000000e+00 & 0.000000e+00 & 0.000000e+00 & 0.000000e+00 & 8.877390e-10 \\ 0.000000e+00 & 0.000000e+00 & 0.000000e+00 & 0.000000e+00 & 0.000000e+00 \\ 0.000000e+00 & 0.000000e+00 & 0.000000e+00 & 0.000000e+00 & 0.000000e+00 \end{pmatrix} \end{equation}\] \[\begin{equation} \myvector{V^T}= \begin{pmatrix} 9.858667e-01 & -6.862581e-02 & 1.131250e-01 & 2.385416e-02 & -9.995558e-02 \\ -8.973915e-02 & -5.978934e-01 & 6.917637e-01 & 1.802947e-01 & 3.513221e-01 \\ 2.199172e-02 & 2.339138e-01 & -1.190638e-01 & 9.531101e-01 & 1.490155e-01 \\ -9.269151e-02 & 5.274638e-01 & 6.849257e-01 & 3.529232e-02 & -4.927692e-01 \\ -1.045870e-01 & -5.521596e-01 & -1.592681e-01 & 2.392914e-01 & -7.755995e-01 \\ \end{pmatrix} \end{equation}\]
この結果からまず\(\myvector{\Lambda}\)が対角行列になっていること、 その対角成分が降順に並んでいることを確かめられる。
From this result, it is confirmed that a diagonal matrix was obtained for \(\myvector{\Lambda}\) whose diagonal components are sorted in a descending order.
次に上記の\(\myvector{U}\), \(\myvector{\Lambda}\), \(\myvector{V}^T\)を用いて \(\myvector{U}\myvector{\Lambda}\myvector{V}^T\)を計算したところ、 以下の結果となった。
Next, \(\myvector{U}\myvector{\Lambda}\myvector{V}^T\) was computed using \(\myvector{U}\), \(\myvector{\Lambda}\), and \(\myvector{V}^T\) obtained above. The result was as follows:
\[\begin{equation} \myvector{U}\myvector{\Lambda}\myvector{V}^T= \begin{pmatrix} -5.542348e-10 & 1.180734e-09 & -5.946389e-10 & -5.832139e-10 & -5.107910e-10 \\ -2.357166e-09 & 1.389110e-09 & 6.061486e-10 & -1.676399e-09 & -8.405894e-10 \\ 9.754172e-10 & -7.404720e-10 & 1.276538e-09 & -9.070130e-10 & -4.890220e-10 \\ 5.082556e-10 & 7.325039e-10 & 6.491890e-10 & 3.288509e-09 & 7.671914e-10 \\ -2.856718e-08 & 1.917195e-09 & -3.393252e-09 & -5.790534e-10 & 2.972070e-09 \\ 1.161411e-09 & 3.587793e-09 & -4.279629e-09 & -4.661626e-10 & -2.133547e-09 \\ -9.078414e-10 & 5.549452e-10 & 7.600898e-10 & 1.573012e-09 & -7.303892e-10 \\ \end{pmatrix} \end{equation}\]
この計算結果と与えた行列\(\myvector{G}\)との差は以下の通りである。
The difference between this result and the given matrix \(\myvector{G}\) was as follows:
\[\begin{equation} \myvector{U}\myvector{\Lambda}\myvector{V}^T-\myvector{G}= \begin{pmatrix} -3.073896e-19 & -2.133711e-18 & 2.324462e-18 & -3.536644e-19 & 9.794455e-19 \\ -1.052869e-18 & -7.214541e-18 & 8.016288e-18 & -1.399970e-19 & 3.607780e-18 \\ 5.176453e-19 & 3.460561e-18 & -3.991493e-18 & -9.186033e-19 & -2.006940e-18 \\ 2.007003e-19 & 1.226150e-18 & -1.614754e-18 & -1.675698e-18 & -1.089717e-18 \\ -1.260381e-17 & -8.398463e-17 & 9.734588e-17 & 2.549924e-17 & 4.960546e-17 \\ 2.759228e-19 & 1.869434e-18 & -2.113169e-18 & -2.060797e-19 & -1.002808e-18 \\ -3.130940e-19 & -2.057880e-18 & 2.434695e-18 & 9.576613e-19 & 1.308821e-18 \\ \end{pmatrix} \end{equation}\]
どの成分で比べても 差\(\myvector{U}\myvector{\Lambda}\myvector{V}^T-\myvector{G}\)は 行列\(\myvector{G}\)そのものに比べて8-9桁ほど小さくなっている。 したがって正しい特異値分解が得られたと判断できる。
For each component, the difference \(\myvector{U}\myvector{\Lambda}\myvector{V}^T-\myvector{G}\) is by 8-9 orders of magnitude smaller than the value of \(\myvector{G}\), indicating that the singular value decomposition given by this function was correct.

◆使用例(Example)