使用 Boost Ublas 时出错 lu_Factorize

Question

我正在尝试制作一个有趣的项目，它制作矩阵幂函数使用 BOOST Ublas.It 与此非常相似 Numpy library power function for matrix 它使用矩阵求幂在对数时间内计算矩阵的 n 次方。计算矩阵的n次方有3种情况-:

如果幂 > 0 使用矩阵求幂直接
如果幂 = 0 检查矩阵是否有逆矩阵（检查使用 lu_factorize）如果是，return 恒等矩阵
如果幂 < 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 倍。你可以观察到但是对于 power <=0，i Get 有时会得到错误的答案。 （我用矩阵和逆矩阵的乘积是单位矩阵的想法检查了这个）

这可能与 lu_factorize 和 lu_substitute 进行某些检查以确保变量类型正确。

由于 lu_factorize 没有可用的文档，我不确定如何使用它。（我刚刚在这里阅读了一个可用的示例 using Ublas lu_factorize and lu_substitute to compute matrix inverse ）。我什至阅读了源代码，但由于缺乏经验，代码掌握不多。

我现在有几个关于我遇到的问题的问题。由于这是我第一次尝试提升，如果我问一些愚蠢的问题，请不要对我太苛刻。所以这些是我的问题 -

错误的答案可能是由于数据类型之间的错误转换或类似的原因。我该如何解决这个问题？我如何确保在每一步都使用正确的数据类型？
当用户输入不正确的类型时，我如何给出错误，我知道我可以使用 boost assert 但那是给出大量难以理解的编译错误。什么是确保输入类型有效的最简单方法。例如我想要如果用户提供输入字符串，则给出错误。你能举个例子吗？
我试过各种方法解决编译错误其中之一是使用 #define BOOST_UBLAS_TYPE_CHECK 0 这有助于绕过数据类型的编译错误，但后来我得到了错误的答案。你能解释一下在这种情况下至少是如何工作的吗？
因为这是我第一次尝试用 boost 做任何事情，我可以理解这肯定可以做得更好。此头文件中必须包含哪些其他内容以确保库标准，例如错误处理，支持多个编译器等？

Answer 1

您的 matrix_power() 函数不 return 来自最终 else 块的值。如果我在该块的末尾添加 return true 那么你的代码运行s 对我来说：

$ 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。

使用 Boost Ublas 时出错 lu_Factorize

Error in using Boost Ublas lu_Factorize

c++

boost

matrix

matrix-multiplication

boost-ublas