GSLを用いて超定方程式群及び行列の奇異値分解(SVD)を解く
18526 ワード
GSLを用いて超定方程式群及び行列の奇異値分解(SVD)を解く
最近,超定方程式群を解く必要がある高動画像(HDR)合成アルゴリズムを学習しているので,GSLを用いてこの問題をどのように解決するかについて,点時間をかけて検討した.
GSLには最小二乗フィッティング(Least-Squares Fitting)に関するアルゴリズムがあり、これらのアルゴリズムの宣言はgsl_fit.hでは、そのままGSLで提供されるgsl_fit_linear関数はこの問題を解決します.しかし、ついでにSVDについてもっと勉強したいと思います.だからgslを直接使っていませんfit_linear関数.
SVD分解のいくつかの基本概念
SVDについては2編の良い科学普及文があります. A Singularly Valuable Decomposition: The SVD of a Matrix We Recommend a Singular Value Decomposition
皆さんに読んでもらうことをお勧めします.この2つの文章はすでに中国語に翻訳されているようです.
SVDとは、マトリックスAmを×nは3つの特殊行列Umに分解する×n 、 Sn×n 、 Vn×nの積.
Am×n=Um×n⋅Sn×n⋅VTn×n
上式のTは行列の転置を表す.分解後のこの3つの行列はまたいくつかの特殊な条件を満たす必要があり,ここでUmは×nとVn×nは直交行列、すなわち、以下を満たす.
UT⋅U=IVT⋅V=I
行列Sは対角行列であり、主対角線上の要素のみが0ではない.
マトリクスUm×n 、 Sn×n 、 Vn×nはいずれも良好な性質を持っているので,このような分解は元の行列Aの性質をよりよく理解するのに役立つ.
例えば、行列Aが満ランク行列である場合、Aは可逆的である.Aの逆は次のように書くことができます.
A−1=V⋅S−1⋅UT
ここでVとUは直交行列であるため、V−1=VT、U−1=UTとなる.Sは対角行列であり,求逆も簡単で,主対角線上の要素ごとに逆数をとるだけである.
GSLにおける相関関数
gslではSVDを計算するためにいくつかの関数が提供されています. gsl_linalg_SV_decompこれは最も基本的なもので、Golub-Reinsch SVDアルゴリズムを使って、一般的にはこれで十分です. gsl_linalg_SV_decomp_modこれは改良されたGolub−Reinsch SVDアルゴリズムであり,M′Nの場合Golub−Reinsch SVDアルゴリズムよりも速い. gsl_linalg_SV_decomp_JAcobiというアルゴリズムはJacobi直交化を用いており,計算結果はGolub−Reinsch SVDアルゴリズムよりも正確であると称されている.
それ以外にもgslがありますlinalg_SV_solve関数.これはSVDの結果を用いて線形代数方程式群を解いたものである.
これらの関数を組み合わせると,線形代数方程式群A
次は関数コードです.
Aが満ランク行列である場合,計算したxは我々の一般的な意味での方程式の解である.
次に、具体的な例を示します.
⎛⎝⎜⎜⎜⎜⎜⎜1.41.63.84.62.62.11.58.08.22.92.11.19.68.40.17.40.75.40.49.99.65.08.88.07.7⎞⎠⎟⎟⎟⎟⎟⎟⋅x=⎛⎝⎜⎜⎜⎜⎜⎜1.11.64.79.10.1⎞⎠⎟⎟⎟⎟⎟⎟
次はテストコードです.
出力結果は次のとおりです.
-5.208566 5.736694 -2.537472 -1.029814 0.968151
1.100000 1.600000 4.700000 9.100000 0.100000
計算結果はまだ正確であることがわかる.
Aの行数が列数より大きい場合に求められるのは最小二乗の意味での解であり,||A⋅x−b||2の最小の解である.次に例を示します.
⎛⎝⎜2314−52⎞⎠⎟⋅x=⎛⎝⎜1136⎞⎠⎟
テストコードは次のとおりです.
計算結果は次のとおりです.
3.090909 1.254545
11.200000 3.000000 5.600000
Aがランクに満たない場合、xは唯一ではない.このとき算出された解の1つ.次に例を示します.
(1224)⋅x=(36)
方程式は簡単で、口算で結果を出すことができます.この方程式の解は:
x=(11)+(−21)⋅t
次に私たちのコードで計算します.
結果は次のとおりです.
-3.400000 3.200000
3.000000 6.000000
(−3.4,3.2)Tは確かに方程式の解であると試算できる.実はSVDで方程式のすべての解を求めることができますが、SとVの値が必要なので、上のlinearSolve_SVD関数は足りません.
次に,SVDに関連する機能をクラスにカプセル化し,SとVの値を抽出するのを容易にする.また,我々の1つのAに複数のxのグループがある場合でも,SVD分解を1回計算するだけで,以下のクラスで多くの計算量を減らすことができる.
ヘッダファイルは次のとおりです.
cppファイルは以下の通りです.
次に、このクラスを使用して、さっきの問題を計算します.
結果は次のとおりです.
S = 5.000000 0.000000
V = -0.447214 -0.894427 -0.894427 0.447214
x = -3.400000 3.200000
Sの2番目の要素は0であることを記し,Vの対応する列(2番目の列)が方程式解の自由ベクトルであることを示した.方程式の解は次のように書くことができます.
x=(−3.43.2)+(−0.8944270.447214)⋅t
この解が正しいことを検証することができます.
また、私が書いたクラスにはtrimVectorS(double abseps)関数も用意されています.この関数を利用して、absepsより小さいすべてのSを直接0に置き換えることができます.この関数が与えられたのは,計算誤差などの影響で,Sの中には0であるはずの項で計算される可能性のある結果が0ではないためである.この関数でこの問題を解決することができます.条件数が大きく,方程式が病的である行列もあり,この関数でもいくつかの問題を解決できる.
はい、先にこんなにたくさん書きます.皆さんに役に立つことを願っています.
最近,超定方程式群を解く必要がある高動画像(HDR)合成アルゴリズムを学習しているので,GSLを用いてこの問題をどのように解決するかについて,点時間をかけて検討した.
GSLには最小二乗フィッティング(Least-Squares Fitting)に関するアルゴリズムがあり、これらのアルゴリズムの宣言はgsl_fit.hでは、そのままGSLで提供されるgsl_fit_linear関数はこの問題を解決します.しかし、ついでにSVDについてもっと勉強したいと思います.だからgslを直接使っていませんfit_linear関数.
SVD分解のいくつかの基本概念
SVDについては2編の良い科学普及文があります.
皆さんに読んでもらうことをお勧めします.この2つの文章はすでに中国語に翻訳されているようです.
SVDとは、マトリックスAmを×nは3つの特殊行列Umに分解する×n 、 Sn×n 、 Vn×nの積.
Am×n=Um×n⋅Sn×n⋅VTn×n
上式のTは行列の転置を表す.分解後のこの3つの行列はまたいくつかの特殊な条件を満たす必要があり,ここでUmは×nとVn×nは直交行列、すなわち、以下を満たす.
UT⋅U=IVT⋅V=I
行列Sは対角行列であり、主対角線上の要素のみが0ではない.
マトリクスUm×n 、 Sn×n 、 Vn×nはいずれも良好な性質を持っているので,このような分解は元の行列Aの性質をよりよく理解するのに役立つ.
例えば、行列Aが満ランク行列である場合、Aは可逆的である.Aの逆は次のように書くことができます.
A−1=V⋅S−1⋅UT
ここでVとUは直交行列であるため、V−1=VT、U−1=UTとなる.Sは対角行列であり,求逆も簡単で,主対角線上の要素ごとに逆数をとるだけである.
GSLにおける相関関数
gslではSVDを計算するためにいくつかの関数が提供されています.
それ以外にもgslがありますlinalg_SV_solve関数.これはSVDの結果を用いて線形代数方程式群を解いたものである.
これらの関数を組み合わせると,線形代数方程式群A
次は関数コードです.
void linearSolve_SVD(const gsl_matrix * A, const gsl_vector * b, gsl_vector * x)
{
int rows = A->size1;
int cols = A->size2;
gsl_vector * work = gsl_vector_alloc (cols);
gsl_vector * S = gsl_vector_alloc (cols);
gsl_matrix * U = gsl_matrix_alloc(rows, cols);;
gsl_matrix * V = gsl_matrix_alloc(cols, cols);
gsl_matrix_memcpy (U, A); // A ， U
gsl_linalg_SV_decomp( U, V, S, work );
gsl_linalg_SV_solve ( U, V, S, b, x );
gsl_vector_free(work);
gsl_vector_free(S);
gsl_matrix_free(V);
gsl_matrix_free(U);
}Aが満ランク行列である場合,計算したxは我々の一般的な意味での方程式の解である.
次に、具体的な例を示します.
⎛⎝⎜⎜⎜⎜⎜⎜1.41.63.84.62.62.11.58.08.22.92.11.19.68.40.17.40.75.40.49.99.65.08.88.07.7⎞⎠⎟⎟⎟⎟⎟⎟⋅x=⎛⎝⎜⎜⎜⎜⎜⎜1.11.64.79.10.1⎞⎠⎟⎟⎟⎟⎟⎟
次はテストコードです.
void test1()
{
double a_data[] = {1.4, 2.1, 2.1, 7.4, 9.6,
1.6, 1.5, 1.1, 0.7, 5.0,
3.8, 8.0, 9.6, 5.4, 8.8,
4.6, 8.2, 8.4, 0.4, 8.0,
2.6, 2.9, 0.1, 9.9, 7.7};
gsl_matrix_view A = gsl_matrix_view_array (a_data, 5, 5);
double b_data[] = {1.1, 1.6, 4.7, 9.1, 0.1};
gsl_vector_view b = gsl_vector_view_array (b_data, 5);
gsl_vector * x = gsl_vector_alloc (5);
linearSolve_SVD(&A.matrix, &b.vector, x);
gsl_vector_fprintf (stdout, x, "%f");
qDebug() << "";
gsl_vector * bb = gsl_vector_alloc (5);
gsl_blas_dgemv (CblasNoTrans, 1, &A.matrix, x, 0, bb);
gsl_vector_fprintf (stdout, bb, "%f");
}出力結果は次のとおりです.
-5.208566 5.736694 -2.537472 -1.029814 0.968151
1.100000 1.600000 4.700000 9.100000 0.100000
計算結果はまだ正確であることがわかる.
Aの行数が列数より大きい場合に求められるのは最小二乗の意味での解であり,||A⋅x−b||2の最小の解である.次に例を示します.
⎛⎝⎜2314−52⎞⎠⎟⋅x=⎛⎝⎜1136⎞⎠⎟
テストコードは次のとおりです.
void test3()
{
double a_data[] = {2, 4,
3, -5,
1, 2};
gsl_matrix_view A = gsl_matrix_view_array (a_data, 3, 2);
double b_data[] = {11, 3, 6};
gsl_vector_view b = gsl_vector_view_array (b_data, 3);
gsl_vector * x = gsl_vector_alloc (2);
linearSolve_SVD(&A.matrix, &b.vector, x);
gsl_vector_fprintf (stdout, x, "%f");
qDebug() << "";
gsl_vector * bb = gsl_vector_alloc (3);
gsl_blas_dgemv (CblasNoTrans, 1, &A.matrix, x, 0, bb);
gsl_vector_fprintf (stdout, bb, "%f");
}計算結果は次のとおりです.
3.090909 1.254545
11.200000 3.000000 5.600000
Aがランクに満たない場合、xは唯一ではない.このとき算出された解の1つ.次に例を示します.
(1224)⋅x=(36)
方程式は簡単で、口算で結果を出すことができます.この方程式の解は:
x=(11)+(−21)⋅t
次に私たちのコードで計算します.
void test4()
{
double a_data[] = {1, 2,
2, 4};
gsl_matrix_view A = gsl_matrix_view_array (a_data, 2, 2);
double b_data[] = {3, 6};
gsl_vector_view b = gsl_vector_view_array (b_data, 2);
gsl_vector * x = gsl_vector_alloc (2);
linearSolve_SVD(&A.matrix, &b.vector, x);
gsl_vector_fprintf (stdout, x, "%f");
qDebug() << "";
gsl_vector * bb = gsl_vector_alloc (2);
gsl_blas_dgemv (CblasNoTrans, 1, &A.matrix, x, 0, bb);
gsl_vector_fprintf (stdout, bb, "%f");
}結果は次のとおりです.
-3.400000 3.200000
3.000000 6.000000
(−3.4,3.2)Tは確かに方程式の解であると試算できる.実はSVDで方程式のすべての解を求めることができますが、SとVの値が必要なので、上のlinearSolve_SVD関数は足りません.
次に,SVDに関連する機能をクラスにカプセル化し,SとVの値を抽出するのを容易にする.また,我々の1つのAに複数のxのグループがある場合でも,SVD分解を1回計算するだけで,以下のクラスで多くの計算量を減らすことができる.
ヘッダファイルは次のとおりです.
#ifndef GSLSINGULARVALUEDECOMPOSITION_H
#define GSLSINGULARVALUEDECOMPOSITION_H
#include cppファイルは以下の通りです.
#include "gsl_SVD.h"
void linearSolve_SVD(const gsl_matrix * A, const gsl_vector * b, gsl_vector * x)
{
int rows = A->size1;
int cols = A->size2;
gsl_vector * work = gsl_vector_alloc (cols);
gsl_vector * S = gsl_vector_alloc (cols);
gsl_matrix * U = gsl_matrix_alloc(rows, cols);;
gsl_matrix * V = gsl_matrix_alloc(cols, cols);
gsl_matrix_memcpy (U, A); // A ， U
gsl_linalg_SV_decomp( U, V, S, work );
gsl_linalg_SV_solve ( U, V, S, b, x );
gsl_vector_free(work);
gsl_vector_free(S);
gsl_matrix_free(V);
gsl_matrix_free(U);
}
int GslSVD::trimVectorS(double abseps)
{
int count = 0;
for(int i = 0; i < S->size; i++)
{
if(fabs(gsl_vector_get(S, i)) < abseps)
{
count ++;
gsl_vector_set(S, i, 0);
}
}
return count;
}
gsl_vector * GslSVD::getVectorS()
{
if(S == NULL) return NULL;
gsl_vector * s = gsl_vector_alloc(S->size);
gsl_vector_memcpy(s, S);
return s;
}
gsl_matrix * GslSVD::getMatrixU()
{
if(U == NULL) return NULL;
gsl_matrix * u = gsl_matrix_alloc(U->size1, U->size2);
gsl_matrix_memcpy(u, U);
return u;
}
gsl_matrix * GslSVD::getMatrixV()
{
if(V == NULL) return NULL;
gsl_matrix * v = gsl_matrix_alloc(V->size1, V->size2);
gsl_matrix_memcpy(v, V);
return v;
}
GslSVD::GslSVD()
{
S = NULL;
U = NULL;
V = NULL;
}
void GslSVD::alloc_suv(int rows, int cols)
{
if( S != NULL )
{
gsl_vector_free(S);
gsl_matrix_free(U);
gsl_matrix_free(V);
}
S = gsl_vector_alloc (cols);
U = gsl_matrix_alloc(rows, cols);
V = gsl_matrix_alloc(cols, cols);
}
int GslSVD::SV_decomp(const gsl_matrix * A)
{
int rows = A->size1;
int cols = A->size2;
gsl_vector * work = gsl_vector_alloc (cols);
alloc_suv(rows, cols);
gsl_matrix_memcpy (U, A); // A ， U
int ret = gsl_linalg_SV_decomp( U, V, S, work );
gsl_vector_free(work);
return ret;
}
int GslSVD::SV_decomp_mod(const gsl_matrix * A)
{
int rows = A->size1;
int cols = A->size2;
gsl_vector * work = gsl_vector_alloc (cols);
gsl_matrix *X = gsl_matrix_alloc(cols, cols);
alloc_suv(rows, cols);
gsl_matrix_memcpy (U, A); // A ， U
int ret = gsl_linalg_SV_decomp_mod( U, X, V, S, work );
gsl_matrix_free(X);
gsl_vector_free(work);
return ret;
}
int GslSVD::SV_decomp_jacobi (gsl_matrix * A)
{
int rows = A->size1;
int cols = A->size2;
alloc_suv(rows, cols);
gsl_matrix_memcpy (U, A); // A ， U
int ret = gsl_linalg_SV_decomp_jacobi( U, V, S );
return ret;
}
int GslSVD::SV_solve(const gsl_vector *b, gsl_vector *x)
{
if(U != NULL)
{
return gsl_linalg_SV_solve (U, V, S, b, x);
}
return -1;
}
GslSVD::~GslSVD()
{
if(S != NULL)
{
gsl_vector_free(S);
gsl_matrix_free(V);
gsl_matrix_free(U);
}
}次に、このクラスを使用して、さっきの問題を計算します.
void test5()
{
double a_data[] = {1, 2,
2, 4};
gsl_matrix_view A = gsl_matrix_view_array (a_data, 2, 2);
GslSVD svd;
svd.SV_decomp(&A.matrix);
puts("S = ");
gsl_vector_fprintf (stdout, svd.getVectorS(), "%f");
puts("
V = ");
gsl_matrix_fprintf (stdout, svd.getMatrixV(), "%f");
double b_data[] = {3, 6};
gsl_vector_view b = gsl_vector_view_array (b_data, 2);
gsl_vector * x = gsl_vector_alloc (2);
svd.SV_solve(b, x);
puts("
x = ");
gsl_vector_fprintf (stdout, x, "%f");
}結果は次のとおりです.
S = 5.000000 0.000000
V = -0.447214 -0.894427 -0.894427 0.447214
x = -3.400000 3.200000
Sの2番目の要素は0であることを記し,Vの対応する列(2番目の列)が方程式解の自由ベクトルであることを示した.方程式の解は次のように書くことができます.
x=(−3.43.2)+(−0.8944270.447214)⋅t
この解が正しいことを検証することができます.
また、私が書いたクラスにはtrimVectorS(double abseps)関数も用意されています.この関数を利用して、absepsより小さいすべてのSを直接0に置き換えることができます.この関数が与えられたのは,計算誤差などの影響で,Sの中には0であるはずの項で計算される可能性のある結果が0ではないためである.この関数でこの問題を解決することができます.条件数が大きく,方程式が病的である行列もあり,この関数でもいくつかの問題を解決できる.
はい、先にこんなにたくさん書きます.皆さんに役に立つことを願っています.