Boost.uBLAS で最小二乗法

boost C++ lapack uBLAS

今回は AX = B の最小二乗解(Least squares solution)，あるいは最小ノルム解(Minimum norm solution)を計算します．lapack を用いてこれを計算するには dgels, dgelss, dgelsd という関数を用いればよく，対応する bindings のヘッダファイルは gels.hpp, gelss.hpp, gelsd.hpp になっています．しかしこれらのヘッダファイルの作者は MSVC 環境でしかテストしていないのか，typename が抜けている箇所がいくつもあります．そのため標準準拠なコンパイラではこれらのヘッダを使用することができません．
この問題はバージョン bindings-20081116 でようやく修正されました．
つまり何が言いたいのかというと，今回のプログラムは bindings-20081116 以降のバージョンを用いないと多くの環境でコンパイルできないということです．注意してください．


AX = B の最小二乗解，あるいは最小ノルム解を計算するには，ムーア・ペンローズの擬似逆行列(Moore-Penrose pseudoinverse)を用いればいいのですが，A x trans(A) あるいは trans(A) x A やその逆行列を直接計算してしまうと余計な誤差(Additional rounding errors)を生じさせてしまいます（そもそも計算する必要すらないわけです）．そこで係数行列 A をQR分解，あるいは特異値分解(SVD)する方法が考えられます．どちらの方法も lapack に実装されており，QR分解して解くのが lapack.dgels, 特異値分解して解くのが lapack.dgelss, lapack.dgelsd になります．lapack.dgelss と lapack.dgelsd の違いは，後者が分割統治法の特異値分解を用いることです．
lapack.dgels は係数行列 A がフルランク(Full rank)であることを想定していることに注意してください．そのため A がランク落ち(Rank deficient)している可能性がある場合は，lapack.dgelss(lapack.dgelsd) を使用してください．

以下に私の実装例を示します．

math.hpp

#ifndef MATH_HPP_20081118
#define MATH_HPP_20081118

#if defined(_MSC_VER) && (_MSC_VER >= 1020)
# pragma once
#endif


#include <boost/numeric/ublas/fwd.hpp>


namespace math {
    struct full_rank {};
    struct rank_deficient {};

    void least_squares(const boost::numeric::ublas::matrix<double>& A, const boost::numeric::ublas::vector<double>& b, boost::numeric::ublas::vector<double>& x, full_rank);
    void least_squares(const boost::numeric::ublas::matrix<double>& A, const boost::numeric::ublas::matrix<double>& B, boost::numeric::ublas::matrix<double>& X, full_rank);
    void least_squares(const boost::numeric::ublas::matrix<double>& A, const boost::numeric::ublas::vector<double>& b, boost::numeric::ublas::vector<double>& x, rank_deficient);
    void least_squares(const boost::numeric::ublas::matrix<double>& A, const boost::numeric::ublas::matrix<double>& B, boost::numeric::ublas::matrix<double>& X, rank_deficient);
}


#endif

math.cpp

#include "math.hpp"
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/vector.hpp>
#define BOOST_NUMERIC_BINDINGS_USE_CLAPACK
#include <boost/numeric/bindings/lapack/gels.hpp>
#include <boost/numeric/bindings/lapack/gelsd.hpp>
#include <boost/numeric/bindings/traits/ublas_matrix.hpp>
#undef  BOOST_NUMERIC_BINDINGS_USE_CLAPACK


void math::least_squares(const boost::numeric::ublas::matrix<double>& A, const boost::numeric::ublas::vector<double>& b, boost::numeric::ublas::vector<double>& x, math::full_rank rank)
{
    BOOST_UBLAS_CHECK(A.size1() == b.size(), boost::numeric::ublas::external_logic());

    boost::numeric::ublas::matrix<double> B(b.size(), 1), X;

    boost::numeric::ublas::column(B, 0).assign(b);
    least_squares(A, B, X, rank);

    x = boost::numeric::ublas::column(X, 0);
}


void math::least_squares(const boost::numeric::ublas::matrix<double>& A, const boost::numeric::ublas::matrix<double>& B, boost::numeric::ublas::matrix<double>& X, math::full_rank)
{
    namespace ublas  = boost::numeric::ublas;
    namespace lapack = boost::numeric::bindings::lapack;

    BOOST_UBLAS_CHECK(A.size1() == B.size1(), ublas::external_logic());

    ublas::matrix<double, ublas::column_major> CA(A), CX((std::max)(A.size1(), A.size2()), B.size2());
    int                                        info;

    ublas::project(CX, ublas::range(0, B.size1()), ublas::range(0, B.size2())).assign(B);

    info = lapack::gels('N', CA, CX, lapack::optimal_workspace());
    BOOST_UBLAS_CHECK(info == 0, ublas::internal_logic());

    X = ublas::project(CX, ublas::range(0, A.size2()), ublas::range(0, B.size2()));
}


void math::least_squares(const boost::numeric::ublas::matrix<double>& A, const boost::numeric::ublas::vector<double>& b, boost::numeric::ublas::vector<double>& x, math::rank_deficient rank)
{
    BOOST_UBLAS_CHECK(A.size1() == b.size(), boost::numeric::ublas::external_logic());

    boost::numeric::ublas::matrix<double> B(b.size(), 1), X;

    boost::numeric::ublas::column(B, 0).assign(b);
    least_squares(A, B, X, rank);

    x = boost::numeric::ublas::column(X, 0);
}


void math::least_squares(const boost::numeric::ublas::matrix<double>& A, const boost::numeric::ublas::matrix<double>& B, boost::numeric::ublas::matrix<double>& X, math::rank_deficient)
{
    namespace ublas  = boost::numeric::ublas;
    namespace lapack = boost::numeric::bindings::lapack;

    BOOST_UBLAS_CHECK(A.size1() == B.size1(), ublas::external_logic());

    ublas::matrix<double, ublas::column_major> CA(A), CX((std::max)(A.size1(), A.size2()), B.size2());
    int                                        info;

    ublas::project(CX, ublas::range(0, B.size1()), ublas::range(0, B.size2())).assign(B);

    lapack::optimal_workspace work;

    info = lapack::gelsd(CA, CX, work);
    BOOST_UBLAS_CHECK(info == 0, ublas::internal_logic());

    X = ublas::project(CX, ublas::range(0, A.size2()), ublas::range(0, B.size2()));
}

math.cpp に AX = B の最小二乗解，あるいは最小ノルム解を計算する関数 least_squares を実装しました．math.cpp をプロジェクトに加え（他のソースコードと一緒にコンパイル），math.hpp を include することにより least_squares が使用可能になります．
least_squares は AX = B の最小二乗解，あるいは最小ノルム解 X を計算します．ただし，A は M x N 行列，X は N x NRHS 行列，B は M x NRHS 行列（NRHS は 1 以上の任意の数）を表します．NRHS が 1 のとき，X と B は列数 1 の行列になりますが，これはベクトルに他なりません（NRHS が 1 のときは，AX = B という式はベクトルを用いて Ax = b になるのです）．uBLAS ではベクトルを表す型もちゃんと用意されているため，引数にベクトルを取る関数(Ax = b を解く関数)も実装しました．コードを見ると分かりますが，この関数は AX = B を解く関数を利用しているだけです．
関数 least_squares は M x N 行列 A を受け取り，M >= N の過剰決定系(Overdetermined system)のときは最小二乗解(Least squares solution)，M < N の不定系(Underdetermined system)のときは最小ノルム解(Minimum norm solution)を計算します．
関数 least_squares の最後の引数に渡す型によって，係数行列 A がフルランク(Full rank)かそうでないか(Rank deficient)を指定します．math::full_rank を渡すと，係数行列 A をフルランクであるとみなし lapack.dgels を用いて計算を行い，math::rank_deficient を渡すと係数行列はランク落ちしている可能性があるとみなし，lapack.dgelsd を用い分割統治法の特異値分解を用い計算を行います．


以下のテストコードで least_squares の動作を確認できます．

test.cpp

#include "math.hpp"
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/matrix.hpp>


#ifdef _MSC_VER
# pragma comment(lib, "libf2c.lib")
# pragma comment(lib, "BLAS.lib")
# pragma comment(lib, "clapack.lib")
#endif


int main()
{
    namespace ublas = boost::numeric::ublas;

    { // 過剰決定系(Overdetermined system)
        ublas::matrix<double> A(6, 2);
        ublas::matrix<double> B(6, 3), X;

        A(0, 0) = 0; A(0, 1) = +2;
        A(1, 0) = 2; A(1, 1) = -1;
        A(2, 0) = 2; A(2, 1) = -1;
        A(3, 0) = 0; A(3, 1) = +1.5;
        A(4, 0) = 2; A(4, 1) = -1;
        A(5, 0) = 2; A(5, 1) = -1;

        B(0, 0) = 1; B(0, 1) = +4; B(0, 2) = 1;
        B(1, 0) = 1; B(1, 1) = +1; B(1, 2) = 2;
        B(2, 0) = 1; B(2, 1) = -1; B(2, 2) = 1;
        B(3, 0) = 1; B(3, 1) = +3; B(3, 2) = 2;
        B(4, 0) = 1; B(4, 1) = +1; B(4, 2) = 1;
        B(5, 0) = 1; B(5, 1) = -1; B(5, 2) = 1;

        { // A がフルランク(Full rank)の場合
            math::least_squares(A, B, X, math::full_rank());

            // 以下結果の表示
            // 計算結果を http://publib.boulder.ibm.com/infocenter/clresctr/vxrx/topic/com.ibm.cluster.essl43.guideref.doc/am501_f10d032.html と比較してみましょう
            std::cout << X << std::endl;
        }

        { // A がランク落ち(Rank deficient)している可能性がある場合
            math::least_squares(A, B, X, math::rank_deficient());

            // 以下結果の表示
            // 計算結果を http://publib.boulder.ibm.com/infocenter/clresctr/vxrx/topic/com.ibm.cluster.essl43.guideref.doc/am501_f10d032.html と比較してみましょう
            std::cout << X << std::endl;
        }
    }

    { // 不定系(Underdetermined system)
        ublas::matrix<double> A(3, 4);
        ublas::matrix<double> B(3, 4), X;

        A(0, 0) = 0.5; A(0, 1) = +0.5; A(0, 2) = 0.5; A(0, 3) = +0.5;
        A(1, 0) = 0.5; A(1, 1) = -1.5; A(1, 2) = 0.5; A(1, 3) = -1.5;
        A(2, 0) = 1.0; A(2, 1) = +1.0; A(2, 2) = 0.0; A(2, 3) = +1.0;

        B(0, 0) = 1; B(0, 1) = +1; B(0, 2) = 1.0; B(0, 3) = 0;
        B(1, 0) = 1; B(1, 1) = -1; B(1, 2) = 2.5; B(1, 3) = 1;
        B(2, 0) = 1; B(2, 1) = +1; B(2, 2) = 3.0; B(2, 3) = 0;

        { // A がフルランク(Full rank)の場合
            math::least_squares(A, B, X, math::full_rank());

            // 以下結果の表示
            // 計算結果を http://publib.boulder.ibm.com/infocenter/clresctr/vxrx/topic/com.ibm.cluster.essl43.guideref.doc/am501_f10d034.html と比較してみましょう
            std::cout << X << std::endl;
        }

        { // A がランク落ち(Rank deficient)している可能性がある場合
            math::least_squares(A, B, X, math::rank_deficient());

            // 以下結果の表示
            // 計算結果を http://publib.boulder.ibm.com/infocenter/clresctr/vxrx/topic/com.ibm.cluster.essl43.guideref.doc/am501_f10d034.html と比較してみましょう
            std::cout << X << std::endl;
        }
    }

    return 0;
}
  • VC9
  • g++ (GCC) 3.4.4 (cygming special)

でのコンパイルを確認しています．
実行してみると，正しく解を計算できていることが分かります．
エラーが出てコンパイルできない場合は，前回の記事を見て正しくインストールできているか確認してみてください．