Matrix — 行列

概要

Matrix<T> は sangi の動的サイズ密行列テンプレートクラスである。

  • 動的サイズ、行優先 (row-major) 格納
  • 全演算子オーバーロード — 行列積、行列-ベクトル積含む
  • 行列式、逆行列、転置、ノルム、トレース
  • ビュー — 部分行列、行、列のコピーなし参照
  • ファクトリ関数 — zeros, identity, 回転行列
  • 静的サイズ版: StaticMatrix<T, Rows, Cols>

コンストラクタ

Matrix();
Matrix(size_type rows, size_type cols);
Matrix(size_type rows, size_type cols, const T& value);
Matrix(std::initializer_list<std::initializer_list<T>> init);
シグネチャ説明
Matrix()空行列 ($0 \times 0$)
Matrix(rows, cols)サイズ指定 (ゼロ初期化)
Matrix(rows, cols, value)サイズ指定 + 全要素を value で初期化
Matrix({{...}, {...}})2D 初期化子リスト
引数型説明
rowssize_type行数
colssize_type列数
valueconst T&初期値
 
Matrix<double> A(3, 3, 0.0);         // 3x3 零行列
Matrix<double> B = {{1, 2}, {3, 4}}; // 2x2

要素アクセス

reference operator()(size_type row, size_type col);
const_reference operator()(size_type row, size_type col) const;
reference at(size_type row, size_type col);
const_reference at(size_type row, size_type col) const;
pointer data() noexcept;
const_pointer data() const noexcept;
メソッド戻り値型説明
A(i, j)T& / const T&Debug 時 assert チェック
A.at(i, j)T& / const T&境界チェック付き (例外送出)
A.data()T* / const T*生ポインタ (行優先連続配置)

サイズ・容量

size_type rows() const noexcept;
size_type cols() const noexcept;
size_type size() const noexcept;
bool empty() const noexcept;
bool is_square() const noexcept;
void resize(size_type rows, size_type cols);
void resize(size_type rows, size_type cols, const value_type& value);
void reshape(size_type rows, size_type cols);
void clear() noexcept;
メソッド戻り値型説明
rows()size_type行数
cols()size_type列数
size()size_type総要素数 ($= \text{rows} \times \text{cols}$)
empty()bool空か
is_square()bool正方行列か
resize(r, c)voidリサイズ
reshape(r, c)void形状変更 (要素数一致が必要)
clear()void空にする

算術演算子

// Matrix + Matrix
Matrix<T> operator+(const Matrix<T>& lhs, const Matrix<T>& rhs);
Matrix<T> operator-(const Matrix<T>& lhs, const Matrix<T>& rhs);
Matrix<T> operator*(const Matrix<T>& lhs, const Matrix<T>& rhs);

// Matrix + Scalar
Matrix<T> operator*(const Matrix<T>& mat, const Scalar& scalar);
Matrix<T> operator*(const Scalar& scalar, const Matrix<T>& mat);
Matrix<T> operator/(const Matrix<T>& mat, const Scalar& scalar);
Matrix<T> operator-(const Matrix<T>& mat);  // 単項マイナス

// Matrix + Vector
Vector<T> operator*(const Matrix<T>& mat, const Vector<T>& vec);

// 複合代入
Matrix& operator+=(const Matrix& rhs);
Matrix& operator-=(const Matrix& rhs);
Matrix& operator*=(const T& scalar);
Matrix& operator/=(const T& scalar);

// operator^ (転置・エルミート転置・逆行列・冪乗)
Matrix operator^(int n) const;
演算戻り値型説明
A + BMatrix<T>行列加算 (同サイズ)。式テンプレートにより一時オブジェクト不要
A - BMatrix<T>行列減算
A * BMatrix<T>行列積 ($A \in \mathbb{R}^{m \times n}$, $B \in \mathbb{R}^{n \times p}$ $\to$ $C \in \mathbb{R}^{m \times p}$)
A * s, s * AMatrix<T>スカラー倍
A * vVector<T>行列-ベクトル積 $y = Ax$
-AMatrix<T>全要素の符号反転
A ^ 'T'Matrix<T>転置 $A^\top$
A ^ 'H'Matrix<T>エルミート転置 $A^*$ (共役転置)
A ^ (-1)Matrix<T>逆行列 $A^{-1}$
A ^ nMatrix<T>冪乗 $A^n$ (二分累乗法)。$n < 0$ は $(A^{-1})^{|n|}$

注意: C++ の ^ 演算子は優先順位が低い。式中では括弧で囲むこと: (A ^ 'T') * B

auto と式テンプレート: A + B や s * A などの演算結果を auto で受けると、 実際の型は Matrix<T> ではなく式テンプレートの中間ノード型になる。 この中間型を logm() や sqrtm() など Matrix<T> を 要求する関数に渡すとコンパイルエラーになる。

// NG: A の型が MatBinOp<...> になり logm() の型推論が失敗
auto A = matmul(B.transpose(), B) + Matrix<double>::identity(n);
auto L = logm(A);  // コンパイルエラー

// OK: Matrix<double> で受けると式が評価される
Matrix<double> A = matmul(B.transpose(), B) + Matrix<double>::identity(n);
auto L = logm(A);  // OK

式の結果を他の関数に渡す場合は Matrix<T> で明示的に受けること。

Matrix<double> A = {{1, 2}, {3, 4}};
Matrix<double> B = {{5, 6}, {7, 8}};
auto C = A * B;  // {{19, 22}, {43, 50}}

Vector<double> v = {1, 2};
auto w = A * v;  // {5, 11}

auto At   = A ^ 'T';   // 転置
auto Ainv = A ^ (-1);   // 逆行列
auto A3   = A ^ 3;      // A*A*A

行列操作

Matrix transpose() const;
T determinant() const;
[[nodiscard]] Matrix inverse() const;
T trace() const;
T norm() const;
void zero();
void fill(const T& value);
void identity();
static Matrix identity(size_type size);
メソッド戻り値型説明
transpose()Matrix転置 $A^\top$
determinant()T行列式 $\det(A)$ (LU 分解, $O(n^3)$)
inverse()Matrix逆行列 $A^{-1}$ (Gauss-Jordan, $O(n^3)$)
trace()Tトレース $\operatorname{tr}(A) = \sum_i a_{ii}$
norm()TFrobenius ノルム $\|A\|_F = \sqrt{\sum_{ij} a_{ij}^2}$
zero()void全要素をゼロに
fill(value)void全要素を value に
identity()void単位行列にする (正方行列のみ)
identity(n)Matrix$n \times n$ 単位行列 (静的)
 
Matrix<double> A = {{1, 2}, {3, 4}};
std::cout << A.determinant() << std::endl;  // -2
auto Ainv = A.inverse();
auto I = A * Ainv;  // ≈ 単位行列

低レベル自由関数

出力先を引数で指定する形式。ComputePolicy テンプレートパラメータでカスタマイズ可能。

template<...> void add(const MatrixA& a, const MatrixB& b, MatrixC& c);
template<...> void subtract(const MatrixA& a, const MatrixB& b, MatrixC& c);
template<...> void scale(const MatrixA& a, const Scalar& alpha, MatrixB& b);
template<...> void multiply(const MatrixA& a, const MatrixB& b, MatrixC& c);
template<...> void multiply_vector(const MatrixA& a, const VectorX& x, VectorY& y);
template<...> void transpose(const MatrixA& a, MatrixB& b);
template<typename MatrixA> auto trace(const MatrixA& a) -> typename MatrixA::value_type;
template<typename MatrixA> auto frobenius_norm(const MatrixA& a) -> typename MatrixA::value_type;
template<typename MatrixA> auto determinant(const MatrixA& a) -> typename MatrixA::value_type;
template<typename T> Matrix<T> inverse(const Matrix<T>& a);

行・列操作・ビュー

template<typename VectorType> void get_row(size_type row_idx, VectorType& dest) const;
template<typename VectorType> void get_column(size_type col_idx, VectorType& dest) const;
template<typename VectorType> void set_row(size_type row_idx, const VectorType& src);
template<typename VectorType> void set_column(size_type col_idx, const VectorType& src);
MatrixView<MatrixBase> submatrix(size_type r, size_type c, size_type nr, size_type nc);
ConstMatrixView<MatrixBase> submatrix(size_type r, size_type c, size_type nr, size_type nc) const;
MatrixView<MatrixBase> row(size_type row_idx);
MatrixView<MatrixBase> column(size_type col_idx);
メソッド戻り値型説明
get_row(i, dest)void$i$ 行目をベクトルにコピー
get_column(j, dest)void$j$ 列目をベクトルにコピー
set_row(i, src)void$i$ 行目を設定
set_column(j, src)void$j$ 列目を設定
submatrix(r, c, nr, nc)MatrixView部分行列ビュー (コピーなし)
row(i)MatrixView$i$ 行目のビュー
column(j)MatrixView$j$ 列目のビュー

ファクトリ関数

template<typename T> Matrix<T> zeros(std::size_t rows, std::size_t cols);
template<typename T> Matrix<T> ones(std::size_t rows, std::size_t cols);
template<typename T> Matrix<T> eye(std::size_t size);
template<typename T> Matrix<T> rotation2D(T angle);
template<typename T> Matrix<T> rotationX(T angle);
template<typename T> Matrix<T> rotationY(T angle);
template<typename T> Matrix<T> rotationZ(T angle);
template<typename T> Matrix<T> hconcat(const Matrix<T>& a, const Matrix<T>& b);
template<typename T> Matrix<T> vconcat(const Matrix<T>& a, const Matrix<T>& b);
template<typename T> [[nodiscard]] Matrix<T> power(const Matrix<T>& base, unsigned int exponent);
関数戻り値型説明
zeros<T>(r, c)Matrix<T>零行列
ones<T>(r, c)Matrix<T>全要素 $1$
eye<T>(n)Matrix<T>$n \times n$ 単位行列
rotation2D(θ)Matrix<T>2D 回転行列 ($2 \times 2$)
rotationX(θ)Matrix<T>3D X 軸回転 ($3 \times 3$)
rotationY(θ)Matrix<T>3D Y 軸回転
rotationZ(θ)Matrix<T>3D Z 軸回転
hconcat(A, B)Matrix<T>水平結合 $[A \mid B]$
vconcat(A, B)Matrix<T>垂直結合 $\begin{bmatrix}A \\ B\end{bmatrix}$
power(A, n)Matrix<T>行列冪 $A^n$ (二分累乗法, $O(n^3 \log k)$)
 
auto I = eye<double>(3);
auto R = rotationZ<double>(M_PI / 4);  // Z 軸 45° 回転

出力

template<typename T> std::ostream& operator<<(std::ostream& os, const Matrix<T>& mat);
template<typename T> std::string to_string(const Matrix<T>& mat, std::size_t precision = 4);
template<typename T> void print_matrix(const Matrix<T>& mat,
    std::ostream& os = std::cout, std::size_t width = 12, std::size_t precision = 6);
関数戻り値型説明
operator<<std::ostream&ストリーム出力
to_string(A, prec)std::string文字列変換
print_matrix(A, os, w, p)void整形出力 (幅・精度指定)

LaTeX 出力

行列やベクトルを LaTeX 形式の文字列に変換するフリー関数。メンバ関数ではなく、sangi 名前空間のフリー関数として提供される。

行列の LaTeX 出力

// style: "pmatrix" (丸括弧), "bmatrix" (角括弧), "vmatrix" (縦線), "matrix" (括弧なし)
template<typename T>
std::string toLatex(const Matrix<T>& mat,
    const std::string& style = "pmatrix", std::size_t precision = 6);

ベクトルの LaTeX 出力 (列ベクトル)

template<typename T>
std::string toLatex(const Vector<T>& v,
    const std::string& style = "pmatrix", std::size_t precision = 6);

パラメータ

引数説明
styleLaTeX 環境名。"pmatrix" = 丸括弧、"bmatrix" = 角括弧、"vmatrix" = 縦線、"matrix" = 括弧なし。デフォルト "pmatrix"
precision小数点以下の桁数。デフォルト 6

使用例

Matrix<double> A({{1, 2}, {3, 4}});
std::cout << toLatex(A) << std::endl;
// \begin{pmatrix}
// 1 & 2 \\
// 3 & 4
// \end{pmatrix}

Vector<double> v({1, 2, 3});
std::cout << toLatex(v, "bmatrix") << std::endl;
// \begin{bmatrix}
// 1 \\
// 2 \\
// 3
// \end{bmatrix}

StaticMatrix<T, Rows, Cols>

静的サイズ版。スタック上に固定サイズで確保。動的版 Matrix<T> と同等の API を提供する。

StaticMatrix<double, 3, 3> A;
constexpr auto r = A.static_rows();  // 3
constexpr auto c = A.static_cols();  // 3

固有メソッド

static constexpr size_type static_rows();
static constexpr size_type static_cols();
StaticMatrix<T, Cols, Rows> transpose() const;
T determinant() const requires (Rows == Cols);
T trace() const requires (Rows == Cols);
static StaticMatrix identity() requires (Rows == Cols);

StaticMatrix 同士の乗算は StaticMatrix<T, M, P> = StaticMatrix<T, M, N> * StaticMatrix<T, N, P> のように コンパイル時にサイズが検証される。

使用例

基本的な行列の作成と演算

#include <math/core/matrix.hpp>
#include <iostream>
using namespace sangi;

int main() {
    // initializer_list construction
    Matrix<double> A = {{1, 2, 3},
                         {4, 5, 6},
                         {7, 8, 10}};

    // transpose, trace, determinant
    auto At = A.transpose();
    std::cout << "trace(A) = " << A.trace() << std::endl;        // 16
    std::cout << "det(A)   = " << A.determinant() << std::endl;  // -3

    // matrix multiplication
    Matrix<double> B = {{1, 0}, {0, 1}, {1, 1}};
    auto C = A * B;  // 3x2 result
    print_matrix(C);

    // scalar operations
    auto D = 2.0 * A - A;  // == A

    // factory functions
    auto I = eye<double>(3);
    auto R = rotationZ<double>(M_PI / 6);  // Z-axis 30 deg rotation
}

連立方程式を solve() で解く

#include <math/core/matrix.hpp>
#include <math/linalg/solve.hpp>
#include <iostream>
using namespace sangi;

int main() {
    // Solve Ax = b
    Matrix<double> A = {{2, 1, -1},
                         {-3, -1, 2},
                         {-2, 1, 2}};
    Vector<double> b = {8, -11, -3};

    auto x = solve(A, b);
    std::cout << "x = " << x << std::endl;  // {2, 3, -1}

    // verify: A * x == b
    auto residual = A * x - b;
    std::cout << "residual norm = " << norm(residual) << std::endl;  // ~0
}

Matrix<Rational> でヒルベルト行列の厳密な逆行列

#include <math/core/matrix.hpp>
#include <math/core/mp/Rational.hpp>
#include <iostream>
using namespace sangi;

int main() {
    const int n = 4;
    Matrix<Rational> H(n, n);
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < n; ++j)
            H(i, j) = Rational(1, i + j + 1);

    // exact inverse -- no floating-point error
    auto Hinv = H.inverse();
    auto I = H * Hinv;

    // verify: every diagonal element is exactly 1
    for (int i = 0; i < n; ++i)
        std::cout << "I(" << i << "," << i << ") = " << I(i, i) << std::endl;
    // Output: 1/1 for all

    print_matrix(Hinv);
}

関連する数学的背景

以下の記事では、行列演算の基盤となる数学的概念を解説している。