【C++】【Eigen】SciPyのspdiags関数（っぽいもの）を実装する

C++,Eigen

先日C++でプログラムを作成していたところ、PythonのSciPyのspdiags()のようなものが必要になったので、Eigenを使って同様の関数を実装してみました。

SciPyのspdiags()関数とは

SciPyのspdiags()関数は、要素を対角に配置した行列を疎行列として返す関数です。dataであらわされる行列の行をdiagsで指定される位置に配置します。

spdiags — SciPy v1.17.0 Manual

少々わかりにくいので実際の例で説明します。dataに対応する行列を

$$
\displaystyle
\begin{pmatrix}
1 & 2 & 3 \\
4 & 5 & 6 \\
7 & 8 & 9
\end{pmatrix}
$$

とします。diagsを$(-1, 0, 1)$とすると、1行目$(1, 2, 3)$は$-1$の位置に配置されます。$-1$とは対角要素に対して、下三角側に1つずらした位置となります。2行目$(5, 6, 7)$は$0$の位置です。$0$の位置は対角要素です。3行目$(7, 8, 9)$は$1$の位置に配置されます。$1$の位置は対角要素に対して、上三角側に1つずらした位置です。位置を行列形式で表すと以下のようになります。

$$
\displaystyle
\begin{pmatrix}
0 & 1 & 2 \\
-1 & 0 & 1 \\
-2 & -1 & 0
\end{pmatrix}
$$

これを基に実際に配置すると次のようになります。

$$
\displaystyle
\begin{pmatrix}
4 & 8 & 0 \\
1 & 5 & 9 \\
0 & 2 & 6
\end{pmatrix}
$$

ちなみに、3行目$(7, 8, 9)$が8から始まっている理由ですが、実際には以下のように配置されるためです。

$$
\displaystyle
\begin{pmatrix}
\color{lightgray}7 & \color{lightgray}0 & \color{lightgray}0 \\
4 & 8 & 0 \\
1 & 5 & 9 \\
0 & 2 & 6 \\
\color{lightgray}0 & \color{lightgray}0 & \color{lightgray}3
\end{pmatrix}
$$

Eigenによる実装

疎行列への要素の配置

疎行列に要素を配置するには、Tripletを使うのが便利です。Tripletは要素の位置と値を持ったオブジェクトで、std::vector<Triplet>に複数保持できます。これをsetFromTriplets()することで、疎行列を一度に作成できます。

insert()を用いて1要素ずつ配置する方法もありますが、要素数が多数となりがちな疎行列においては、Tripletを使う方法のほうが高速なようです。

実装

Eigenを使い、C++でこれと似たような関数を実装します。

ソースコードはGitHubにも置いてあります。Eigenと同じようにヘッダファイルのインクルードのみで動作します。

GitHub上には、diagsがstd::vector<int>のバージョンと、std::initializer_list<int>のバージョンも用意しました。

#include <algorithm>
#include <stdexcept>
#include <vector>
#include <Eigen/Eigen>

// SciPyのspdiags()関数っぽいもの
//
// data : 対角に配置する要素を列に持つ行列
// diags: dataの各列について配置する位置を指定します。
//        (0: 対角、正の値: 上三角側、負の値：下三角側)
// m    : 出力する疎行列の行数（省略時または負の値はdataに一致）
// n    : 出力する疎行列の列数（省略時または負の値はdataに一致）
template <typename Derived>
Eigen::SparseMatrix<typename Derived::PlainObject::Scalar>
Spdiags(const Eigen::MatrixBase<Derived> & data, const Eigen::VectorXi & diags, const int m = -1, const int n = -1)
{
    using Scalar = typename Derived::PlainObject::Scalar;
    using T = Eigen::Triplet<Scalar>;

    if(diags.size() > data.rows()) {
        throw std::invalid_argument("diags size is larger than rows of data.");
    }

    int row_size, column_size;
    Eigen::SparseMatrix<Scalar> a;
    std::vector<T> triplets;

    row_size = (m <= 0) ? (int)data.rows() : m;
    column_size = (n <= 0) ? (int)data.cols() : n;

    a.resize(row_size, column_size);
    triplets.resize(row_size * diags.size());

    for (int k = 0; k < diags.size(); k++) {
        int start_index = std::max(0, diags(k));
        int end_index = std::min((int)data.cols(), column_size);

        for(int i = start_index; i < end_index; i++) {
            if(i - diags(k) < row_size && i < column_size) {
                triplets.push_back(T(i - diags(k), i, data(k, i)));
            }
        }
    }

    a.setFromTriplets(triplets.begin(), triplets.end());

    return a;
}

使い方

#include <iostream>
#include "spdiags.h"

int main()
{
    Eigen::MatrixXd a(3, 3);
    Eigen::VectorXi d1(3);

    a << 1, 2, 3,
         4, 5, 6,
         7, 8, 9;
    d1 << -1, 0, 1;

    std::cout << Spdiags(a, d1) << std::endl;

    Eigen::MatrixXd b(3, 4);
    Eigen::VectorXi d2(3);

    b << 1,  2,  3,  4,
         5,  6,  7,  8,
         9, 10, 11, 12;
    d2 << 0, 1, 2;

    std::cout << Spdiags(b, d2, 9, 9) << std::endl;

    return 0;
}

出力結果

Nonzero entries:
(4,0) (1,1) (8,0) (5,1) (2,2) (9,1) (6,2)

Outer pointers:
0 2 5  $

4 8 0
1 5 9
0 2 6

Nonzero entries:
(1,0) (6,0) (2,1) (11,0) (7,1) (3,2) (12,1) (8,2) (4,3)

Outer pointers:
0 1 3 6 9 9 9 9 9  $

1 6 11 0 0 0 0 0 0
0 2 7 12 0 0 0 0 0
0 0 3 8 0 0 0 0 0 
0 0 0 4 0 0 0 0 0
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0
  • Pythonの出力結果
>>> from scipy.sparse import spdiags
>>> a = [[1, 2, 3], [4, 5, 6], [7, 8, 9]]
>>> d1 = [-1, 0, 1]
>>> spdiags(a, d1).toarray()
array([[4, 8, 0],
       [1, 5, 9],
       [0, 2, 6]])
>>> b = [[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]]
>>> d2 = [0, 1, 2]
>>> spdiags(b, d2, 9, 9).toarray()
array([[ 1,  6, 11,  0,  0,  0,  0,  0,  0],
       [ 0,  2,  7, 12,  0,  0,  0,  0,  0],
       [ 0,  0,  3,  8,  0,  0,  0,  0,  0],
       [ 0,  0,  0,  4,  0,  0,  0,  0,  0],
       [ 0,  0,  0,  0,  0,  0,  0,  0,  0],
       [ 0,  0,  0,  0,  0,  0,  0,  0,  0],
       [ 0,  0,  0,  0,  0,  0,  0,  0,  0],
       [ 0,  0,  0,  0,  0,  0,  0,  0,  0],
       [ 0,  0,  0,  0,  0,  0,  0,  0,  0]])

Pythonの結果と出力が一致していることがわかります。

動作を確認した環境

  • OS
    1. Windows 11
      1. Microsoft C/C++ Optimizing Compiler 19.44.35213 for x64
      2. MinGW-W64 GCC version 13.2.0
      3. Clang version 20.1.0
    2. Linux
      1. GCC version 12.2.0
  • Eigen version 3.4.0

まとめ

Eigenを使用してSciPyのspdiags()関数のようなものを自作してみました。入力行列dataの行方向の要素をベクトルdiagsで指定した位置に対角方向に配置し、疎行列を作成して返します。

ちなみにMATLABにもspdiags()関数がありますが、Pythonが行方向の要素を対角に配置するのに対して、列方向の要素を対角に配置します。MATLAB風にするには、ソースコードを変える必要がありますのでご注意ください。

C++,Eigen

Posted by izadori