リアルタイム制御向け C++行列ライブラリ

Claude A

2024年11月4日 15:49

目的

行列計算をマイコンなどに実機実装したい

カルマンフィルタ、適応制御、モデル予測制御

しかし、世にあるライブラリは、

ランタイムエラーを回避できない（ゼロ割、マイナス平方根）
解の精度を調整できないため、メモリ消費、計算時間が大きい

という問題があります。
そこで、リアルタイム制御に適した行列計算ライブラリを作ることにしました。
完成品は以下にあります。

コンセプト

密行列、対角行列、スパース行列を使い分けられる
反復計算の回数調整が可能
ゼロ割、マイナス平方根などを回避する機能を実装

以前に作成しました PythonからC++へ自動変換するツールで、Numpyの行列計算コードから簡単に置き換えられるようにインターフェースを作り込みました。

説明動画

詳細については以下の動画で説明していますのでご覧ください。

スパース行列の定義方法を変更（2024/12/01）

スパース行列の定義方法を変更しました。詳しくはこちらの動画をご確認ください。

行列の定義方法について（2025/01/27）

C++テンプレートクラスの性質上、行列の定義が何だかんだ面倒くさいです。
そこで、もっと簡単に定義できる関数を作りました。

auto A = make_DenseMatrix<3, 3>(1.0, 2.0, 3.0, 5.0, 4.0, 6.0, 9.0, 8.0, 7.0);

auto B = make_DiagMatrix<3>(1.0, 2.0, 3.0);

auto C =
    make_SparseMatrix<SparseAvailable<ColumnAvailable<true, false, false>,
                                      ColumnAvailable<true, false, true>,
                                      ColumnAvailable<false, true, true>>>(
        1.0, 3.0, 8.0, 2.0, 4.0);

このように、「make_～～」と書いて、最低限のテンプレートと引数を記述するだけで定義できるようにしました。

密行列の場合はテンプレート引数に行数、列数を指定し、関数引数に要素を1行1列、1行2列、…、M行N列の順番に入れます。
この時、関数引数のデータ型で行列要素のデータ型が決まるので注意してください。

対角行列の場合は、行数 = 列数なので、テンプレート引数に指定するのはMだけです。
関数引数には対角行列の要素を順番に入れます。

スパース行列は、テンプレート引数に"SparseAvailable"を書いて入れます。ここは手間ですが省略できません。
関数引数には0でない要素を順番に入れます。

スパース行列の定義について

スパース行列はCompressed Sparse Row（CSR）方式を用いています。（別名、Compressed Row Storage、CRS形式とも言います）

一つの行列を三つの配列で表現します。

value = {1, 3, 8, 2, 4} // 値
row_index = {0, 0, 2, 1, 2} // 列インデックス
row_pointer = {0, 1, 3, 5} // 行Headポインタ

この三つの配列で表現される元の行列は、以下のようになります。

このスパース行列を用いて、密行列と「スパース行列 * 密行列」の計算を行うと、以下のように0である要素の乗算、加算を行わないで効率的に計算できます。

Matrix<T, M, K> Y;

for (std::size_t i = 0; i < K; i++) {
  for (std::size_t j = 0; j < M; j++) {
    T sum = static_cast<T>(0);
    for (std::size_t k = A.row_pointer(j); k < A.row_pointer(j + 1); k++) {
      sum += A.value(k) * B(A.row_index(k), i);
    }
    Y(j, i) = sum;
  }
}

複素数定義方法について改善（2024/11/04）

二つ目の動画の最後で、複素数の定義がbase_matrix側にあるため、使おうとすると、

#include "python_numpy.hpp"

using namespace PythonNumpy;

auto eigen_values_matrix =
      Matrix<DefDiag, Base::Matrix::Complex<double>, 3>(eigen_values_c.matrix);

のように、「Base::Matrix::Complex」と書かなければなりませんでした。上位レイヤーのPythonNumpyでは、下位レイヤーのBase::Matrixを使わないで書けるようにしたかったのですが、複素数だけはできない、というか、やろうとすると面倒でできずにいました。

いろいろ調べて分かりましたが、以下のように名前空間の置き換え機能を使えば、自然に下位レイヤーの機能を上位レイヤーでラッパーできることに気付きました。

#include "base_matrix.hpp"

namespace PythonNumpy {

template <typename T> using Complex = Base::Matrix::Complex<T>;

} // namespace PythonNumpy

これによって、以下のようにBase::Matrixを書かずに使えるようになりました。

#include "python_numpy.hpp"

using namespace PythonNumpy;

auto eigen_values_matrix =
      Matrix<DefDiag, Complex<double>, 3>(eigen_values_c.matrix);

行列の要素にアクセスする方法（2024/12/02）

本ライブラリでは、行列の要素にアクセスする方法は二つあります。

() を使ってアクセス

行列の変数名にかっこを付けて、行数、列数を指定する方法です。

【密行列の場合】

/* 定義 */
Matrix<DefDense, double, 3, 3> A({ { 1, 2, 3 }, {5, 4, 6}, {9, 8, 7} });

/* アクセス */
double a_value = A(1, 2); // get
A(1, 2) = a_value;        // set

【対角行列の場合】

/* 定義 */
Matrix<DefDiag, T, 3> B({ 1, 2, 3 });

/* アクセス */
double b_value = B(1); // get
B(1) = b_value;        // set

対角行列では、対角の要素を配列として格納していますので、上記のように1を指定すると、2行2列目の要素にアクセスできます。

【スパース行列の場合】

/* 定義 */
Matrix<DefSparse, T, 3, 3,
    SparseAvailable<
    ColumnAvailable<true, false, false>,
    ColumnAvailable<true, false, true>,
    ColumnAvailable<false, true, true>>
    > C({ 1, 3, 8, 2, 4 });

/* アクセス */
double c_value = C(2); // get
C(2) = c_value;        // set

スパース行列では、要素の配列、ここでは{ 1, 3, 8, 2, 4 }に直接アクセスします。上記の例では、2を指定すると、8の値にアクセスできます。

get, setを使ってアクセス

get, setでは、三つの行列タイプ全てで行数、列数を指定します。