我正在尝试制作一个有趣的项目,它制作一个矩阵幂函数使用BOOST Ublas。它与Numpy库的矩阵幂函数非常相似。它使用矩阵幂运算来计算对数时间内矩阵的n次幂。它有三种情况来计算矩阵的n次方-:
- 如果幂>0,则直接使用矩阵幂
- 如果幂=0,请检查矩阵是否有反矩阵(请检查Using lu_factorize),如果是,则返回单位矩阵
- 如果功率<0找到逆(如果存在)然后用矩阵求幂
我擅长算法和实现,但这是当我第一次使用任何开源库时,我想要学习这一点,以便我最终能够为boost做出贡献。
这是我的头文件
// 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 matrixn");
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 inversen");
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 matrixn");
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 << "nnn";}
函数在幂>0时工作得很好,因为它使用矩阵直接求幂。它的工作速度比重复乘法快得多,我已经看到了大约1000次迭代,这和使用循环的运行时间相差100-1000倍。你可以观察到但为了权力<0,我有时会得到错误的答案(我使用矩阵和逆的乘积是单位矩阵的想法来检查这一点)
这可能与lu_factorize和lu_substitute,后者具有某些检查以确保变量类型正确。
由于他们没有关于lu_factorize的文档,我不知道如何使用它。
我现在有几个关于我正在经历的问题的问题。由于这是我第一次尝试助推,如果我问一些愚蠢的问题,请不要对我太苛刻。这些是我的问题-
- 错误的答案可能是由于数据类型之间的转换不正确,或者类似的原因。我该如何解决此问题?如何确保每一步都使用正确的数据类型
- 当用户输入不正确的类型时,我该如何给出错误?我知道我可以使用boost断言,但这是给出了大量无法理解的编译错误。有什么确保输入类型有效的最简单方法。例如,我想要如果用户给出一个字符串进行输入,则会给出一个错误。你能举个例子吗
我尝试了各种方法来解决编译错误其中之一是使用#define BOOST _UBLAS_TYPE_CHECK 0这有助于绕过数据类型的编译错误,但后来我得到了错误的答案。你能解释一下至少在这种情况下是如何工作的吗?
由于这是我第一次尝试从助推中获得任何东西,我可以理解这肯定可以做得更好。这个头文件中还必须包含哪些内容才能确保库标准错误处理,支持多个编译器等?
matrix_power()
函数不会从最后一个else块返回值。如果我在该块的末尾添加return true
,那么您的代码将为我运行:
$ clang++ --version
Apple LLVM version 7.0.0 (clang-700.1.76)
Target: x86_64-apple-darwin14.5.0
Thread model: posix
$ clang++ -std=c++11 -Wall -Wextra -g -I/opt/local/include lu.cpp -o lu
$ ./lu
[3,3]((-0.0644531,-0.00195312,0.861328),(1.09896,-0.0677083,-11.1406),(-0.9375,0.0625,9.4375))
我还能够删除BOOST_UBLAS_TYPE_CHECK
定义,而不会出现编译时(或运行时)错误。我使用的是Boost 1.59。