🍴

fortran呼び出しのlapackで線形最小二乗法を実装する

2024/03/02に公開

C++

lapack

lstsq

tech

c言語はfortranの関数を呼び出すことができます。ここではlapackeなどを使わず、直接lapackを呼び出して、線形最小二乗法を実現するコードの例を示します。

おそらく多くの場合、この記事で紹介している実装方法はあまりに面倒なためアンチパターンといえるでしょう。著者はCのインターフェースであるlapacke.hや、boostのbinding、Eigenなどのライブラリを利用されることをおすすめいたします。実際に線形最小二乗法をboostのublasとbindingを用いたて実装したコードが以下にあります。今回紹介するコードの設計はこのコードを参考に作成しました。
https://yano.hatenadiary.jp/entry/20081118/1227001319

しかしさまざまな事情により、生のlapack呼び出しが必要になることがあるかもしれません。この記事はそんなときのための参考資料にしていただければと思います。

線形最小二乗法とは

$m\times n$ の行列 $A$ と $m$ 列のベクトル $b$ があるとき、 $x$ を $m$ 列のベクトルとして、以下の線形方程式が成り立ちます。

Ax = b

この方程式について、行列 $A$ の列が行よりも多い $m>n$ が成り立つ場合の問題を線形最小二乗法と呼びます。 $m=n$ が成り立つとき、すなわち連立方程式の問題のとき、 $A$ の逆行列を $b$ にかけたものが解になりますが、 $m>n$ の場合についても同様の演算子があり、疑似逆行列 $A^{\dagger}$ を $b$ に作用させた値が最小二乗

\min || Ax - b ||

を満たす解となります。詳しくは疑似逆行列のWikipediaが詳しいでしょう。

したがって線形最小二乗法は疑似逆行列を求める問題に帰着され、行列の問題として処理することができるようになります。実際のルーチンでは直接疑似逆行列を求めない方法を用いていますが、大まかなコンセプトは上記のアイデアです。

lapackでの実装

lapackにおいて、線形最小二乗法を計算するルーチンは、いくつかあります。

?GELSS: 特異値分解を使って、優決定系または劣決定系の連立1次方程式Ax=Bに対する最小ノルム最小二乗解を求める
?GELSD: Householder変換で元の問題を二重対角最小二乗問題にして、それを解いたうえでさらにHouseholder変換で戻してもとの最小二乗解を求める
?GELSY: 完全直交分解を使って最小二乗解を求める

ここで最初の?は数値精度に関するオプションを表していて、私の知る限りのlapackでは、

S: 32bit 浮動小数点数
D: 64bit 浮動小数点数
C: 64bit 複素数浮動小数点数
Z: 128bit 複素数浮動小数点数

が実装されています。

ではどのルーチンを使うべきかという疑問がわくでしょう。しっかりやるならご自身の問題に応じて、速度を調べるのがよろしいかと思いますが、ここでは参考として、上記の3つのルーチンのpythonのラッパーであるscipy.linalg.lstsqの説明を紹介します。

Which LAPACK driver is used to solve the least-squares problem. Options are 'gelsd', 'gelsy', 'gelss'. Default ('gelsd') is a good choice. However, 'gelsy' can be slightly faster on many problems. 'gelss' was used historically. It is generally slow but uses less memory.

線形最小二乗法を解くのにどのlapackドライバーを使うとよいだろうか？オプションは、gelsd、gelsy、gelssがある。デフォルト（gelsd）は良い選択肢である。しかし、gelsyは多くの問題で若干速いに違いない。gelssは古くから使われてきた。これは多くのケースで遅いが、メモリ利用量が少なくすむ。（筆者訳）

まとめると困ったらgelsdを利用するとよいそうです。そこでscipyの記述に従い、本記事ではdgelsdをc++言語からfortranを呼び出し、利用する例を紹介します。

環境構築

lapackのうまい環境構築について別にまとめた記事を作成しました。そちらを参照していただけたらと思います。以下の環境構築では、lapackeをinstallしていますので、本記事と合わせてご参照いただければと思います。

dgelsdのapi:

$Ax=b$ の最小二乗法について、 $A$ が $m\times n\quad ( m>n)$ の行列で、 $b$ が $m$ 行のベクトルとします。
dgelsdは以下の引数を持ちます。

int *m: $A$ の行数
int *n: $A$ の列数
int *nrhs:解きたい $b$ の数、今回は1を指定します。
double *A:行列 $A$ （列優先の配列で渡す）
double *lda: $A$ のleading dimension size(mを代入すればok)
double *b:ベクトル $b$
double *ldb: $b$ のleading dimension size(mを代入すればok)
double *s:結果の特異値が格納される
double *rcond:精度-1を指定するとマシンイプシロンで求めてくれる
double *rank:Aの有効なrank
double *work:作業領域
int *lwork:作業領域
int *iwork:作業領域
int *info:成功かどうかのフラグ

ここで作業領域に何を与えればよいのだろうかと疑問を抱くかもしれません。APIリファレンスを読むと行列サイズなどから指定される式のサイズを与えればよいそうですが、すこし自分で計算式を書くのは面倒です。実は便利なことには、lworkに-1を格納してdgelsdルーチンを実行するとこの作業領域を求めて、work変数、iwork変数に格納して返してくれるという機能があります！

よって作業手順は次のようになります：

列優先で配列A,b,sを作成
lwork=-1を指定してdgelsdルーチンを呼ぶ
その結果を受けて、workスペースの配列を生成
dgelsdルーチンを再び呼び、計算
結果を出力

fortranのコードをc++で呼ぶ際の注意点

基本的にfortranもcもc++も共有ライブラリを言語に関係なく共有できるため、fotranから生成された共有ファイルをcから何も問題なく呼び出し、コンパイルすることは可能なのですがc++から呼び出す場合は注意する必要があります。それは、C++はコンパイルする際に名前マングリングと呼ばれる処理が行われることです。つまりc++では低レイヤーにおける関数名にもともとつけた名前にプレフィックスが付与された名前が付けるようになっています。その関係で共有ファイルの関数について、単に前方宣言を書いて呼び出すことができません。それを回避するにはextern "C" {}で囲ってあげる必要があります。これについては以下の記事が詳しいでしょう：

実装

まずは最小二乗法を実行するヘッダーファイルです：

least_square.hpp

#pragma once 
#include <vector>
#include <exception>
#include <iostream>

extern "C" {
  int dgelsd_(int *m, int *n, int *nrhs, double *A, int *lda, double *b, int *ldb, double *s, double *rcond, int *rank, double *work, int *lwork, int *iwork, int *info);
}

// A,b (A: m x n where m > n,b: m) -> min|| Ax -y ||^2 -> x, singluar value
namespace LeastSquare {
  struct result {
    std::vector<double> xs;
    std::vector<double> singular_values;
  };
  // require: A.size() == m A[:].size() == n  b.size() == m
  // caution: Calling this function, b would be changed.
  result calc(std::vector<std::vector<double>> &A, std::vector<double> &b);
};

LeastSquare::result LeastSquare::calc(std::vector<std::vector<double>> &A, std::vector<double> &b)
{
  if (A.size() != b.size())
  {
    throw std::runtime_error("The row of A must be same as the the row of b.");
  }
  int m = A.size(), n = A[0].size(), nrhs = 1, lda = m, ldb = std::max(m,n), info, rank;
  double rcond = -1.0;
  std::vector<double> singular_values(std::min(m,n));
  std::vector<double> A_col_major(m*n);
  
  for (int i=0; i < m; ++i)
  {
    for (int j=0; j < n; ++j)
    {
      A_col_major[j*m +i] = A[i][j];
    }
  }

  // check workspace:
  int lwork=-1, iwork; double work_size;
  dgelsd_(&m, &n, &nrhs, A_col_major.data(), &lda, b.data(), &ldb, singular_values.data(), &rcond, &rank, &work_size, &lwork, &iwork, &info);
  // actual run:
  lwork = static_cast<int>(work_size);
  double *work = new double[lwork];
  int *iwork2 = new int[iwork];
  dgelsd_(&m, &n, &nrhs, A_col_major.data(), &lda, b.data(), &ldb, singular_values.data(), &rcond, &rank, work, &lwork, iwork2, &info);
  delete[] work;
  delete[] iwork2;
  // make_result:
  if (info != 0)
  {
    std::runtime_error("calculation failed at dgelsd_");
  }
  std::vector<double> xs(n);
  std::copy(b.begin(), b.begin() + n, xs.data());
  return {xs, singular_values};
}

続いて、このヘッダーファイルをテストするコードです。

check.cpp

#include "leqst_square.hpp"
#include <vector>
#include <random>
#include <iostream>
int main()
{
  std::mt19937 rng(30);
  std::uniform_real_distribution<double> rand(0, 1);
  int m = 4, n = 2;
  std::vector<std::vector<double>> A(m, std::vector<double>(n));
  std::vector<double> b(m);
  for (int i=0; i < A.size(); ++i)
  {
    for (int j=0; j < A[0].size(); ++j)
    {
      A[i][j] = rand(rng);
    }
    b[i] = rand(rng);
  }

  std::cout << "A:\n";
  for (auto row: A)
  {
    for (auto val: row)
    {
      std::cout << val << ",";
    }
    std::cout << "\n";
  }
  std::cout << "b:\n";
  for (auto val: b)
  {
    std::cout << val << ",";
  }
  LeastSquare::result result = LeastSquare::calc(A, b);
  std::cout << "\nsingular value:\n";
  for (auto val : result.singular_values)
  {
    std::cout << val << ",";
  }
  std::cout << std::endl;
  std::cout <<"x:\n";
  for (auto val: result.xs)
  {
    std::cout << val << ",";
  }
  std::cout << "\n";
  return 0;
}

コンパイルは次のように行います。lapackの環境構築編で紹介した方法でlapackをインストールされている場合は、

g++ check.cpp -o check -llapack

でできますが、そうではない場合は、

 find / -name "liblapack.so*" 2>/dev/null

で検索しましょう。きっとシステムに1つや2つはlapackが入っているに違いありません！！（入っていない場合は前述の環境構築編などを参考にしていただき、installしていただければと思います。
私の場合の一例として、

/usr/lib/x86_64-linux-gnu/liblapack.so

がありました。その場合、

g++ check.cpp /usr/lib/x86_64-linux-gnu/liblapack.so -o check

でコンパイルできるでしょう。もちろんこの共有ファイルが入っているディレクトリに`LD_LIBRARY_PATHや/etc/ld.so.confが指定されていることが実行には必須です。

さて実行すると、

A:
0.916978,0.218257,
0.643062,0.706414,
0.927394,0.696651,
0.440224,0.177017,
b:
0.409313,0.600071,0.0235988,0.249151,
singular value:
1.78835,0.416345,
x:
0.28574,0.169201,

と表示されるでしょう。（seedが30としています。）

これはpythonで簡単に確認することができます。scipy.linalg.lstsqにこの行列とベクトルを入れてxと特異値が出てくるかを調べましょう：

check.py

import numpy as np
from scipy.linalg import lstsq

A = np.array([
    0.916978, 0.218257,
    0.643062, 0.706414,
    0.927394, 0.696651,
    0.440224, 0.177017,
]).reshape(4, 2)
b = np.array([
    0.409313, 0.600071, 0.0235988, 0.249151,
])

result = lstsq(A, b, cond=-1)
print("singular value:", result[3])
print("x:", result[0])

実行すると、

$python check.py
singular value: [1.78835275 0.41634495]
x: [0.28574033 0.16920022]

となるでしょう。これで検証も完了です。

終わりに

ここまで読んでいただき、ありがとうございました。今回はlapackのdgelsdをcppからlapackeなどのラッパーライブラリを用いずに直接呼び出し、最小二乗法を実行するコードを紹介しました。lapackは非常に有名なライブラリであるのに対して、こうしたドキュメントが少なくとも日本語ではインターネットには転がっていなかったため、本記事を作成しました。ですがよく考えるとこんなニッチなことをしている人がいないということなのかもしれません。

だいたいlapackeを使えよって話ですよね。ちなみに私の場合は数値計算コードを書いていて、後からlapackを利用したいとなったのですが、元の数値計算コードの依存関係のライブラリが利用していたlapackがlapackeなしでシステムにインストールされていた関係でlapackeが使えませんでした。いやきっとうまくやれば使えるのかもしれませんが、とにかく今回はlapackを生で使うことにしました。

皆さまの研究・開発が一層進むことを願いまして、本記事を終えようと思います。