BOOST Ublas を使用して行列べき乗関数を作成する楽しいプロジェクトを作成しようとしています。このNumpy ライブラリの行列のべき乗関数と非常によく似ています。行列指数を使用して、行列の n 乗を対数時間で計算します。行列の n 乗を計算するには、次の 3 つのケースがあります。
- べき乗 > 0 の場合は、行列の累乗を直接使用します
- べき乗 = 0 の場合、行列に逆行列があるかどうかをチェックし (lu_factorize を使用してチェック)、ある場合は単位行列を返します。
- べき乗< 0 の場合、逆行列 (存在する場合)を見つけてから、 行列累乗を使用します。
私はアルゴリズムと実装に長けていますが、オープン ソース ライブラリを使用するのはこれが初めてです。これを学びたいので、最終的にブーストに貢献できます。
これは私のヘッダーファイルです
// Distributed under the Boost Software License, Version 1.0. (See
// accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
//
#ifndef BOOST_UBLAS_POW_HPP
#define BOOST_UBLAS_POW_HPP
#include <iostream>
#include <vector>
#include <iomanip>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/triangular.hpp>
#include <boost/multiprecision/cpp_int.hpp>
#include <boost/mpl/set.hpp>
#include <boost/mpl/assert.hpp>
#include <boost/multiprecision/cpp_int.hpp>
using namespace boost::numeric::ublas;
namespace boost { namespace numeric { namespace ublas {
typedef permutation_matrix<std::size_t> pmatrix;
template< typename T,typename T2 >
matrix<T> matrix_power(const matrix<T> input, T2 exponent)
{
matrix<T> resultant=input;
BOOST_ASSERT_MSG(input.size1()==input.size2(),"Not a square matrix\n");
if(exponent>0)
resultant=matrix_exponent(input,exponent);// this where you could directly use matrix exponentiation
else if(exponent==0)
{
pmatrix pm(input.size1());
matrix<T> A(input);
BOOST_ASSERT_MSG(lu_factorize(A, pm)==0,"Attempted to compute a power of 0 in a matrix without inverse\n");
resultant.assign(identity_matrix<T> (input.size1()));// if the matrix is invertible output identity matrix for the 0t hpower
}
else {
matrix<T> A(input);
pmatrix pm(A.size1());
BOOST_ASSERT_MSG(lu_factorize(A, pm)==0,"Attempted to compute inverse in a singular matrix\n");
resultant.assign(identity_matrix<T> (A.size1()));
lu_substitute(A, pm, resultant);
resultant=matrix_exponent(resultant,std::abs(exponent));
}// in last case first we compute the inverse of the matrix and then we do matrix exponentiation
return resultant;
}
template< typename T,typename T2 >
matrix<T> matrix_exponent(matrix<T> base,T2 exponent)
{
matrix<T> resultant=base;
exponent-=2;
while(exponent>0)
{
if(exponent%2==1)resultant=prod(resultant,base);
exponent=exponent >> 1;
base=prod(base,base);
}
return resultant;
}
}
}
}
#endif
これを使用してこのヘッダーファイルをテストしています
#include "power.hpp"
using namespace boost::numeric::ublas;
using namespace std;
typedef permutation_matrix<std::size_t> pmatrix;
int main()
{
matrix<double> m (3, 3);
for(int i=0;i<3;i++)for(int j=0;j<3;j++)m(i,j)=3*i+j+8;
m(0,0)=11;
matrix<double> c(3, 3);
int h=matrix_power(m,c,-1);// c stores -1 power of m
if(h)// h tells whether the power exists or not
std:: cout << c << "\n\n\n";}
関数は、累乗 > 0 に対して完全に正常に機能します。これは、行列の累乗を直接使用するためです。反復乗算よりもはるかに高速に動作します。約1000回の反復で、これの実行時間とループの使用時間が100〜1000倍異なることがわかりました。あなたはそれを観察することができますが、べき乗 <=0,i の場合、時々間違った答えが得られます。（行列と逆行列の積が単位行列であるという考えを使用してこれを確認しました）
これはおそらく、 変数の型が正しいことを確認する特定のチェックを行うlu_factorize および lu_substitute と関係があります。
lu_factorize のドキュメントがないため、使用方法がわかりません。( Ublas の lu_factorize と lu_substitute を使用して逆行列を計算する例をここで読んだだけです)。ソース コードも読んでいますが、コードは経験不足のためあまり理解していません。
現在発生している問題についていくつか質問があります。これはブーストの私の最初の試みなので、私が愚かなことを尋ねても、私に厳しくしないでください. これらは私の質問です -
- 間違った答えは、おそらくデータ型の間の変換が正しくないか、または同様のものによるものです。どうすればこれを解決できますか? すべてのステップで正しいデータ型を確実に使用するにはどうすればよいですか?
- ユーザーが間違った型を入力したときにエラーを出すにはどうすればよいですか。ブーストアサートを使用できることはわかっていますが、理解できないコンパイルエラーが大量に発生しています。入力タイプが有効であることを確認する最も簡単な方法は何ですか。たとえば、ユーザーが入力用の文字列を指定した場合にエラーを発生させたいとします。これについて例を挙げていただけますか？
コンパイル エラーを解決するためにさまざまな方法を試しましたが、そのうちの 1 つは#define BOOST_UBLAS_TYPE_CHECK 0を使用することでし た。この場合、少なくともそれがどのように機能するか説明できますか？
ブーストから何かを作るのはこれが初めての試みなので、これは確かにもっと良くできることは理解できます. エラー処理、複数コンパイラのサポートなどのライブラリ標準を保証するために、このヘッダー ファイルに存在する必要がある他のものは何ですか?