C++で行列計算実装

はじめに

C++初心者が練習も兼ねて行列計算をC++使って自前で実装してみました。この実装で行列積や逆行列、行列式計算等基本的な行列計算はカバーできています。簡単に拡張できるため特殊な行列演算をしたい場合などに使えます。

なお今回は勉強のために行列計算を自分で実装してみましたが、C++にはEigenという行列計算ライブラリがあるらしいので早く済ませるならそちらを使うのがよいと思います。
Eigen使い方

実現すること

簡単のためN行xM列の形の行列に形式を制限して実装します。ベクトルを扱いたい場合はNx1の行列として定義します。

下記のような計算を実現します。

//A, B, CをNxM行列とする
//和
C = A + B; //要素ごとの和
C = A + 3; //Aの各要素に3足す
C = 3 + A; 
C = A - B; //要素ごとの差
C = A - 3; //Aの各要素から3引く
C = 3 - A; 

//AをNxM行列, BをMxS行列, CをNxS行列とする
//積
C = A * B; //行列積

//A, CをNxN行列とする
C = A * 3; //Aの各要素を３倍
C = 3 * A;
C = -A; //Aの各要素の符号反転

//除算
C = A / 3; //Aの各要素を3で割る

//転置
C = A.T();

//逆行列
C = A.inv();

//行列式
C = A.det();

//Aの中身を表示
A.show(); 
//もしくは
show(A);

Maクラスの宣言

Maというクラスをこの行列計算クラスとします。中身はvector<vector<double>>です。したがってこれから実装するもろもろの関数の宣言を省くと下記になります。

class Ma{
  private:
    vector<vector<double>> vec; //このオブジェクトが保持する行列
  public:
    Ma(vector<vector<double>> x){vec=x;}; //変換コンストラクタ
    operator vector<vector<double>>(){return vec;}; //変換関数
}

ここでvecがこのオブジェクトが扱う行列で、コンストラクタで2次元vector型からMa型に変換できるようにしておきます。また、逆にMa型からvectorに戻すためにoperator vector<vector<double>>()でvecを返すようにしておきます。

AをMa型、aをvector<vector<double>>型とするとこれによりA=(Ma)aといった変換やa=(vector<vector<double>>)Aといった変換が可能になりました。

以降は各演算の説明に入ります。

和算の実装

足し算を実装できることがこのプログラムの基本になります。これができれば文法的には他の演算もすべて実装できるようになります。

行列同士の和

足し算を実装するには、+演算子をオーバーロードするので、宣言は下記のようにoperator()+という関数として追加します。

class Ma{
  private:
    vector<vector<double>> vec;
  public:
    Ma(vector<vector<double>> x){vec=x;};
    operator vector<vector<double>>(){return vec;};
    Ma operator+(Ma x); //追加
}

関数の中身自体はとても簡単です。オブジェクトが保持する行列をA、+演算子の引数として入力された行列をBとします。また結果を保持する行列をCとします。Aの行ごと及び列ごとで２重のfor文を回してCにAとBの各要素ごとの和を代入、最後にCを返すだけです。関数の実装は下記です。なお、本来AとBの行列形式が一致していることを確認するエラー処理があるべきですが一旦省略します。最後にソースコードを貼り付けているのでそちらではエラー処理も入れました。

Ma Ma::operator+(Ma x){
  vector<vector<double>> A=vec;
  vector<vector<double>> B=x.vec;
  vector<vector<double>> C;
  C = A;

  for(int i=0; i < A.size(); i++){
    for(int j=0; j < A[i].size(); j++){
      C[i][j] = A[i][j] + B[i][j];
    }
  }
  Ma y(C);
  return y;
}

行列とスカラー値の和

次に行列の各要素にスカラー値を加算する関数を実装します。

A+3のように左オペランドが行列だった場合は行列同士の和と同じ関数で引数の型をdouble等に変更するだけで実現できますが、3+Aのように左オペランドが数値の場合そうはいきません。
そこでこれをfriend関数として定義します。friend関数はクラスのメンバ関数ではありませんが、クラスのprivateメンバにアクセスできる特別な関数です。

class Ma{
  private:
    vector<vector<double>> vec;
  public:
    Ma(vector<vector<double>> x){vec=x;};
    operator vector<vector<double>>(){return vec;};
    Ma operator+(Ma x); 
    friend Ma operator+(Ma X, double y); //追加 （これはメンバ関数としても定義可能）
    friend Ma operator+(double y, Ma X); //追加
}

関数の実装は行列同士の和計算よりさらに簡単です。さきほどと同様に各要素を取り出してから各要素にスカラー値を加算して結果を返すだけです。

Ma operator+(Ma X, double y){
  vector<vector<double>> x = X.vec;
  vector<vector<double>> A(x.size(), vector<double>(x[0].size(),y));
  Ma Y(A);
  return X + Y;
}
Ma operator+(double y, Ma X){
  return X + y;
}

行列同士の差

行列同士の差は和の実装時における+を-に読み替えるだけなので詳細の説明は省略します。関数の中身は下記になりました。

Ma Ma::operator-(Ma x){
  vector<vector<double>> A=vec;
  vector<vector<double>> B=x.vec;
  vector<vector<double>> C;
  C = A;

  for(int i=0; i < A.size(); i++){
    for(int j=0; j < A[i].size(); j++){
      C[i][j] = A[i][j] - B[i][j];
    }
  }
  Ma y(C);
  return y;
}

行列とスカラー値の差

和と同様にfriend関数で定義・実装します。

Ma operator-(Ma X, double y){
  vector<vector<double>> x = X.vec;
  vector<vector<double>> A(x.size(), vector<double>(x[0].size(),y));
  Ma Y(A);
  return X - Y;
}
Ma operator-(double y, Ma X){
  return -X + y;
}

ここで行列Xの符号を反転する計算-Xが出てきていますが、これはまだ定義していませんでした。下記で定義します。

Ma Ma::operator-(){
  vector<vector<double>> A=vec;
  vector<vector<double>> C = A;
  for(int i=0;i<A.size();i++){
    for(int j=0;j<A[0].size();j++){
      C[i][j] = -A[i][j];
    }
  }
  this->vec = C;
  return *this;
}

前置演算子-をオーバーロードしています。前置なのに引数なしの形式で書く理由は謎ですがこの書き方でうまくいきました。

積計算の実装

行列積は実装がすこしややこしいので先に数値と行列の積を実装します。

数値と行列の積

行列の各要素に特定の数値をかける処理は和の計算のときとほとんど同じです。下記のようにfriend関数で実装します。

Ma operator*(Ma X, double y){
  vector<vector<double>> x = X.vec;
  vector<vector<double>> A(x.size(), vector<double>(x[0].size(),y));
  for(int i=0;i<x.size();i++){
    for(int j=0;j<x[0].size();j++){
      A[i][j] = x[i][j] * y;
    }
  }
  return (Ma)A;
}
Ma operator*(double y, Ma X){
  return X * y;
}

数値による行列の除算

除算の場合も同様です。ただしメンバ関数として実装しました。なお、先に実装した積を利用してもっと簡潔に書く方法もあると思います。

Ma Ma::operator/(double x){
  vector<vector<double>> A=vec;
  for(int i=0;i < A.size();i++){
    for(int j=0;j<A[0].size();j++){
      A[i][j] = vec[i][j] / x;
    }
  }
  return A;
}

行列積

行列積を計算します。NxM行列とMxS行列の積はNxS行列になるはずで、第１オペランドの行列をA, 第２オペランドの行列をBとするとAのi行目とBのj列目の内積をとって結果を結果格納行列のi,j成分に格納します。

計算方法の詳細は解説しているサイトが山ほどあるので下記など他サイトに譲ります。
行列積の計算方法

上記で説明したことをコードにすると下記のようになりました。

Ma Ma::operator*(Ma x){
  vector<vector<double>> A=vec;
  vector<vector<double>> B=x.vec;
  vector<vector<double>> C(A.size(), vector<double>(B[0].size(),0));

  for(int i=0; i < A.size(); i++){
    for(int s=0;s<B[0].size();s++){
      C[i][s] = dot(A[i],get_column(B,s));
    }
  }
  Ma y(C);
  return y;
}

なお、上記のために内積計算する関数double dot(vector<double> a,vector<double> b)と行列のs番目の列を取り出す関数vector<double> get_dolumn(vector<vector<double>> B, int s)を定義しました。

vector<double> Ma::get_column(vector<vector<double>> x,int n){
  vector<double> y(x.size());
  for(int i=0;i<x.size();i++){
    y[i] = x[i][n];
  }
  return y;
}
double Ma::dot(vector<double> x, vector<double> y){
  double z = 0;
  for(int i=0;i<x.size();i++){
    z += x[i] * y[i];
  }
  return z;
}

なお、今回の実装では*で行列積を計算することにしました。*で同じサイズの行列の要素ごとの積（アダマール積）を計算したほうが嬉しい場合もあります。今回のクラスを拡張してアダマール積を導入する場合は別の演算子で実装するのがいいと思います。

転置

転置はNxM行列のオブジェクトに対してMxNの結果行列を作成し、オブジェクトのi,j成分を結果のj,i成分に代入するだけです。実装は下記です。

Ma Ma::T(){
  vector<vector<double>> A=this->vec;
  vector<vector<double>> C(A[0].size(), vector<double>(A.size()));
  for(int i=0;i<A.size();i++){
    for(int j=0;j<A[0].size();j++){
      C[j][i] = A[i][j];
    }
  }
  Ma y(C);
  return y;
}

逆行列計算

逆行列は掃き出し法を使います。これも計算方法の詳細は他に解説サイトがたくさんあるのでここでは省略します。
掃き出し法計算方法

Ma Ma::inv(){
  vector<vector<double>> A(vec.size(), vector<double>(vec[0].size()*2,0));//A[N][N*2]の行列を作る
  int n = vec.size();

  for(int i=0;i<n;i++){
    for(int j=0;j<n;j++){
      A[i][j] = vec[i][j]; //Aの前半列にはオブジェクトの行列を格納
      A[i][j+n] = (i==j)? 1:0; //Aの後半列には単位行列を格納
    }
  }

  for(int i=0;i<n;i++){
    A[i] = div(A[i],A[i][i]); //Aのi行目についてA[i][i]の値で割る
    for(int j=0;j<n;j++){
      if(i==j)
        continue;
      double t = A[j][i];
      for(int k=0;k<n*2;k++){
        A[j][k] = A[j][k] - A[i][k] * t; //Aのj行目についてA[j][i]=0となるようにAのi行目の定数倍を引く
      }
    }
  }
  //この時点でAの前半は単位行列になっている
  vector<vector<double>> B(n, vector<double>(n));
  for(int i=0;i<n;i++){
    for(int j=0;j<n;j++){
      B[i][j] = A[i][j+n]; //Aの後半列を結果行列に格納
    }
  }
  return (Ma)B;
}

ただし、途中でベクトルの各要素を数値で除算する関数divを使いました。

vector<double> Ma::div(vector<double> x, double y){
  for(int i=0;i<x.size();i++){
    x[i] = x[i] / y;
  }
  return x;
}

行列式の計算

行列式は余因子展開する方法で実装しました。余因子展開の詳細は他のサイトを参照してください。
余因子展開計算方法

double Ma::det(){
  vector<vector<double>> A(vec.size(), vector<double>(vec[0].size()));
  A = vec;
  int n=A.size();
  if(n==1)
    return A[0][0];
  else if(n==2)
    return A[0][0] * A[1][1] - A[0][1] * A[1][0]; //要素数２までは簡単なので直接計算

  //0行目で余因子展開
  int sum = 0;
  for(int i=0;i<n;i++){ //A[0][i]で余因子展開する
    vector<vector<double>> B(A.size(), vector<double>(A[0].size()));
    B = A;
    for(int j=0;j<n;j++){
      B[j].erase(B[j].begin()+i);//B[j][i]を消す
    }
    B.erase(B.begin()); //0行目は最後に消す
    sum += A[0][i] * pow(-1,i+2) * ((Ma)B).det(); //Bが余因子行列になっている
  }
  return sum;
}

ただし、累乗を計算する関数powを利用しました。

double Ma::pow(double x, int n){
  int r = 1;
  for(int i=0;i<n;i++){
    r *= x;
  }
  return r;
}

その他補助機能

最後に、行列の中身を簡単に確認するために行列内容を表示する関数や行列のサイズを返す関数を定義します。

行列内容の表示

行列Aの中身を表示するには、A.show()のような書き方やMa::show(A)のような書き方ができるのが便利です。これらを下記で実装しました。

void Ma::show(){
  vector<vector<double>> A=this->vec;
  cout << "\n";
  for(int i=0; i < A.size(); i++){
    for(int j=0; j < A[i].size(); j++){
      cout << A[i][j] << " ";
    }
    cout << "\n";
  }
}

void Ma::show(Ma x){
  x.show();
}

例えば(Ma){{1,2},{3,4}}をMa::show(A)で表示した場合下記のように表示されます。

1 2
3 4

行列の要素数の確認

行列Aの要素数を整数で返す関数sizeを下記で実装しました。A.size()のように呼びます。

int Ma::size(){
  return vec.size() * vec[0].size();
}

行列の形式の確認

NxM行列のNとMが知りたい場合に使います。

vector<int> Ma::shape(){
  vector<int> a = {(int)vec.size(),(int)vec[0].size()};
  return a;
}

おわりに

以上で実装方法の解説は全てです。C++初心者のため稚拙・冗長・高計算量なコードになってしまった部分もあるかもしれません。その点は申し訳ありません。

ソースコード

関数を宣言するヘッダファイルと関数の定義で別々にしてあります。

mat.h

#include <iostream>
#include <vector>
using namespace std;

class Ma{
  private:
    vector<double> get_column(vector<vector<double>> x, int n);
    double dot(vector<double> x, vector<double> y);
    vector<vector<double>> vec;
    vector<double> div(vector<double> x, double y);
    double pow(double x, int n);
  public:
    Ma(vector<vector<double>> x){vec=x;};//変換コンストラクタ
    void show();
    static void show(Ma x);
    Ma operator+(Ma x);
    friend Ma operator+(Ma X, double y);
    friend Ma operator+(double y, Ma X);
    Ma operator*(Ma x);
    friend Ma operator*(Ma X, double y);
    friend Ma operator*(double y, Ma X);
    Ma operator-(Ma x);
    friend Ma operator-(Ma X, double y);
    friend Ma operator-(double y, Ma X);
    Ma operator/(double x);
    Ma operator-();
    Ma T();
    operator vector<vector<double>>(){return vec;};//変換関数
    double operator()(int x, int y){return vec[x][y];};
    int size();
    vector<int> shape();
    Ma inv(); //逆行列
    double det(); //行列式
};

mat.cpp

#include <iostream>
#include <vector>
#include "mat.h"
using namespace std;

double Ma::pow(double x, int n){
  int r = 1;
  for(int i=0;i<n;i++){
    r *= x;
  }
  return r;
}
int Ma::size(){
  return vec.size() * vec[0].size();
}
vector<int> Ma::shape(){
  vector<int> a = {(int)vec.size(),(int)vec[0].size()};
  return a;
}
Ma Ma::operator-(Ma x){
  vector<vector<double>> A=vec;
  vector<vector<double>> B=x.vec;
  vector<vector<double>> C;
  C = A;

  try{
    if(A.size()!=B.size() || A[0].size()!=B[0].size())
      throw "Invalid matrix size! in Ma Ma::operator-(Ma x)";
  }catch(const char* err){
    cout << "Error:" << err << "\n";
    exit(1);
  }

  for(int i=0; i < A.size(); i++){
    for(int j=0; j < A[i].size(); j++){
      C[i][j] = A[i][j] - B[i][j];
    }
  }
  Ma y(C);
  return y;
}
Ma Ma::operator+(Ma x){
  vector<vector<double>> A=vec;
  vector<vector<double>> B=x.vec;
  vector<vector<double>> C;
  C = A;
  
  try{
    if(A.size()!=B.size() || A[0].size()!=B[0].size())
      throw "Invalid matrix size! in Ma Ma::operator+(Ma x)";
  }catch(const char* err){
    cout << "Error:" << err << "\n";
    exit(1);
  }

  for(int i=0; i < A.size(); i++){
    for(int j=0; j < A[i].size(); j++){
      C[i][j] = A[i][j] + B[i][j];
    }
  }
  Ma y(C);
  return y;
}
Ma Ma::operator*(Ma x){
  vector<vector<double>> A=vec;
  vector<vector<double>> B=x.vec;
  vector<vector<double>> C(A.size(), vector<double>(B[0].size(),0));
  try{
    if(A[0].size()!=B.size()){
      cout << "The number of columns in first operand:" << A[0].size() << "\n";
      cout << "The number of rows in second operand:" << B.size() << "\n";
      throw "Invalid matrix size! in Ma Ma::operator*(Ma x)";
    }
  }catch(const char* err){
    cout << "Error:" << err << "\n";
    exit(1);
  }

  for(int i=0; i < A.size(); i++){
    for(int s=0;s<B[0].size();s++){
      C[i][s] = dot(A[i],get_column(B,s));
    }
  }
  Ma y(C);
  return y;
}
Ma Ma::operator/(double x){
  vector<vector<double>> A=vec;
  for(int i=0;i < A.size();i++){
    for(int j=0;j<A[0].size();j++){
      A[i][j] = vec[i][j] / x;
    }
  }
  return A;
}
vector<double> Ma::get_column(vector<vector<double>> x,int n){
  vector<double> y(x.size());
  for(int i=0;i<x.size();i++){
    y[i] = x[i][n];
  }
  return y;
}
double Ma::dot(vector<double> x, vector<double> y){
  double z = 0;
  for(int i=0;i<x.size();i++){
    z += x[i] * y[i];
  }
  return z;
}
void Ma::show(){
  vector<vector<double>> A=this->vec;
  cout << "\n";
  for(int i=0; i < A.size(); i++){
    for(int j=0; j < A[i].size(); j++){
      cout << A[i][j] << " ";
    }
    cout << "\n";
  }
}

void Ma::show(Ma x){
  x.show();
}

Ma Ma::T(){
  vector<vector<double>> A=this->vec;
  vector<vector<double>> C(A[0].size(), vector<double>(A.size()));
  for(int i=0;i<A.size();i++){
    for(int j=0;j<A[0].size();j++){
      C[j][i] = A[i][j];
    }
  }
  Ma y(C);
  return y;
}
vector<double> Ma::div(vector<double> x, double y){
  for(int i=0;i<x.size();i++){
    x[i] = x[i] / y;
  }
  return x;
}
Ma Ma::inv(){
  vector<vector<double>> A(vec.size(), vector<double>(vec[0].size()*2,0));
  try{
    if(vec.size()!=vec[0].size())
      throw "Invalid matrix size. inv() can only be called by square matrix!";
  }catch(const char* err){
    cout << "Error:" << err << "\n";
    exit(1);
  }
  int n = vec.size();
  for(int i=0;i<n;i++){
    for(int j=0;j<n;j++){
      A[i][j] = vec[i][j];
      A[i][j+n] = (i==j)? 1:0;
    }
  }

  for(int i=0;i<n;i++){
    A[i] = div(A[i],A[i][i]);
    for(int j=0;j<n;j++){
      if(i==j)
        continue;
      double t = A[j][i];
      for(int k=0;k<n*2;k++){
        A[j][k] = A[j][k] - A[i][k] * t;
      }
    }
  }
  
  vector<vector<double>> B(n, vector<double>(n));
  for(int i=0;i<n;i++){
    for(int j=0;j<n;j++){
      B[i][j] = A[i][j+n];
    }
  }
  return (Ma)B;
}

double Ma::det(){
  vector<vector<double>> A(vec.size(), vector<double>(vec[0].size()));
  A = vec;
  int n=A.size();
  if(n==1)
    return A[0][0];
  else if(n==2)
    return A[0][0] * A[1][1] - A[0][1] * A[1][0];
  try{
    if(vec.size()!=vec[0].size())
      throw "Invalid matrix size. det() can only be called by square matrix!";
  }catch(const char* err){
    cout << "Error:" << err << "\n";
    exit(1);
  }
  //0行目で余因子展開する
  int sum = 0;
  for(int i=0;i<n;i++){
    vector<vector<double>> B(A.size(), vector<double>(A[0].size()));
    B = A;
    for(int j=0;j<n;j++){
      B[j].erase(B[j].begin()+i);
    }
    B.erase(B.begin()); //1行目は最後に消す
    sum += A[0][i] * pow(-1,i+2) * ((Ma)B).det();
  }
  return sum;
}


Ma Ma::operator-(){//前置-
  vector<vector<double>> A=vec;
  vector<vector<double>> C = A;
  for(int i=0;i<A.size();i++){
    for(int j=0;j<A[0].size();j++){
      C[i][j] = -A[i][j];
    }
  }
  this->vec = C;
  return *this;
}
Ma operator+(Ma X, double y){
  vector<vector<double>> x = X.vec;
  vector<vector<double>> A(x.size(), vector<double>(x[0].size(),y));
  Ma Y(A);
  return X + Y;
}
Ma operator+(double y, Ma X){
  return X + y;
}
Ma operator*(Ma X, double y){
  vector<vector<double>> x = X.vec;
  vector<vector<double>> A(x.size(), vector<double>(x[0].size(),y));
  for(int i=0;i<x.size();i++){
    for(int j=0;j<x[0].size();j++){
      A[i][j] = x[i][j] * y;
    }
  }
  return (Ma)A;
}
Ma operator*(double y, Ma X){
  return X * y;
}
Ma operator-(Ma X, double y){
  vector<vector<double>> x = X.vec;
  vector<vector<double>> A(x.size(), vector<double>(x[0].size(),y));
  Ma Y(A);
  return X - Y;
}
Ma operator-(double y, Ma X){
  return -X + y;
}