1.はじめに
今回はタイトルにもあるように行列演算ライブラリをC++で作りたいと思いました。ライブラリもある中なぜわざわざこのようなことを考えたのかというとただの自己満です。あとは単純に勉強になると思ったからです。また１からといいながら生成AIを使いまくっています。
少しまじめな話をしますと行列の演算ができるようになったらそのライブラリを使っていろいろ遊んでやろうと思っています。といっても今のところはニューラルネットワークを組んで学習させるといったことしか考えていませんが…
2.まずはベクトル
行列を計算しようにもいきなり行列を作成するわけではありません。大体はVectorを使うことが多いと思います(Vectorを利用して行列を計算する人もあまりいない気はしますが)。しかし今回はベクトルの計算から自分で実装していきます。今回はコードを実行するときの効率や厳密なメモリ管理といったところはあまり考えずにひたすら完成を目指します。
3.ベクトルの初期化
ベクトルはテンプレートでサイズを決めてそこから動的にメモリを確保するという方式を採用しました。そこでまずはテンプレート引数でベクトルの型とサイズを受け取ってnewでメモリを確保します。
以下にコンストラクタを示しておきます。
Vector(){
vector = new T[Dim];
}
Vector(T value){
vector = new T[Dim];
Init_vec<T, Dim>::init(vector, value);
}
Vector(T min, T max){
vector = new T[Dim];
Init_vecRan<T, Dim>::init(vector, min, max);
}
Vector(std::initializer_list<T> initList){
if (initList.size() != Dim) {
throw std::runtime_error("Initializer list size does not match vectorsize.");
}
vector = new T[Dim];
std::copy(initList.begin(), initList.end(), vector);
}
初期化の際に{1, 2, 3}のような値を一つずつ直に指定できるようにしたかったのでinitializer_listを使っています。
ここで初期化も行うわけですが今回はテンプレート再帰という方法で初期化を行っています。ここはシンプルにfor文を回すという選択肢もあったかとは思いますが勉強のためでもあるため新しいことに挑戦しました。
ここでテンプレート再帰に関して少し説明します。これは関数による再帰のような感じで構造体の中でテンプレート引数を1ずつ減らしながらテンプレートに代入して自身を呼び出すというものです。そしてとある値の時に終了するというものです。言葉では分かりにくいと思うのでコードを載せておきます。
//定数で初期化するためのテンプレート再帰
template<typename Type, size_t Dim>
struct Init_vec{
static void init(Type *vector, Type value){
vector[Dim - 1] = value;
Init_vec<Type, Dim - 1>::init(vector, value);
}
};
template<typename Type>
struct Init_vec<Type, 0>{
static void init(Type *vector, Type value){}
};
また乱数でも初期化できるようにしました。これはニューラルネットワークを構築するときに重みが必要であることを見越して組み込みました。
//乱数で初期化するためのテンプレート再帰
std::random_device rd;
std::mt19937 gen(rd());
template<typename Type, size_t Dim>
struct Init_vecRan {
static void init(Type *vector, Type min, Type max) {
if constexpr (std::is_integral<Type>::value) {
std::uniform_int_distribution<Type> dis(min, max);
vector[Dim - 1] = dis(gen);
Init_vecRan<Type, Dim - 1>::init(vector, min, max);
} else {
std::uniform_real_distribution<Type> dis(min, max);
vector[Dim - 1] = dis(gen);
Init_vecRan<Type, Dim - 1>::init(vector, min, max);
}
}
};
template<typename Type>
struct Init_vecRan<Type, 0> {
static void init(Type *vector, Type min, Type max) {}
};
これに加えてベクトルにはベクトルの要素へのアクセスと代入そして全体の表示が必要なのでそちらを追加しておきます。要素のアクセスはVectorがconstの場合とそうでない場合の二通りを用意します。そして代入に関しては問題ないかと思いますが表示に関してはベクトルの数字の型が分からないのでstd::coutによる表示を行っております。
T& operator[](size_t i){
if(i >= Dim) throw std::out_of_range("Index out of range.");
return vector[i];
}
const T& operator[](size_t i)const{
if(i >= Dim) throw std::out_of_range("Index out of range.");
return vector[i];
}
Vector& operator=(const Vector& A){
for(size_t i = 0; i < Dim; i++)
this->vector[i] = A.vector[i];
return *this;
}
void print() const {
for (size_t i = 0; i < Dim; ++i) {
std::cout << vector[i] << " ";
}
std::cout << std::endl;
}
4.ベクトルから行列へ
ベクトルから行列を作成するには結論から言うとベクトルのひとつひとつの要素にベクトルを代入すればよいのです。例えば三段ボックスを３つ横にずらっと並べると3x3行列のような感じになると思うのですがそのボックスがここでいうベクトルに当たります。ようするにベクトルを並べれば1次元だったものが2次元として扱えるということです。早速その宣言を見てみましょう。
//行列の確保
Vector<Vector<Type, Col>, Row> matrix;
Vectorは型と大きさを宣言するのですが型の部分がintやdoubleではなくVectorになっています。こうすることでVectorの中身をVectorにすることができました。またRowは行でColは列の大きさを表しています。
以下はコンストラクタなのですがfor文を使ってVectorのそれぞれの要素に入っているVectorを順に初期化していきます。またベクトルの時と同じようにinitializeer_listで{{1, 2}, {3, 4}}のような初期化を可能としています。
Matrix(Type value){
for(int i = 0; i < Row; i++){
matrix[i] = Vector<Type, Col>(value);
}
}
Matrix(Type min, Type max){
for(int i = 0; i < Row; i++){
matrix[i] = Vector<Type, Col>(min, max);
}
}
Matrix(std::initializer_list<std::initializer_list<Type>> initList) {
if (initList.size() != Row) {
throw std::runtime_error("Initializer list row size does not match matrix row size.");
}
int i = 0;
for (auto rowList : initList) {
if (rowList.size() != Col) {
throw std::runtime_error("Initializer list column size does not match matrix column size.");
}
matrix[i++] = Vector<Type, Col>(rowList);
}
}
以上が行列の初期化の原理です。
5.メンバたち
行列を作成することはできました。ですが計算ができなければ作る意味はありません。ということでここからは各演算を実装していきます。とその前にこちらも要素へアクセスできるようにします。アクセスの仕方は２通りあります。１つは普通の２次元配列と同じように[2][3]とかくパターンです。もう１つは(2, 3)とかくパターンです。
Vector<Type, Col>& operator[](size_t i){
return matrix[i];
}
const Vector<Type, Col>& operator[](size_t i)const{
return matrix[i];
}
Type& operator()(size_t i, size_t j){
return matrix[i][j];
}
const Type& operator()(size_t i, size_t j)const{
return matrix[i][j];
}
それではここから演算子を定義していきます。というのも演算子にはもともとこのように使ってくださいと決められているのでMatrix + Matrixのように使ったところでエラーが出るだけです。なのでMatrix用に演算子を定義するわけですね。具体的には先ほどからちらほらと目にしたであろうoperatorです。これは演算子のオーバーロードを行うための命令です。ちなみにオーバーロードとは関数でいうと全く同じ名前の関数でも引数が違えば別のものとして扱われるというものです。A()という関数とA(int x)は別物ということですね。これを演算子でもやってしまおうというのがoperatorなのです。通常+という演算子はintやdoubleといった型を受け取って計算しているわけですがVectorを受け取るタイプの+も作るといった具合です。実際に見たほうが早いかもしれませんね。
Matrix& operator=(const Matrix& A){
for(int i = 0; i < Row; i++){
for(int j = 0; j < Col; j++){
(*this)(i, j) = A(i, j);
}
}
return *this;
}
template<typename Expr>
Matrix& operator=(const Expr& e){
for(int i = 0; i < Row; i++){
for(int j = 0; j < Col; j++){
(*this)(i, j) = e(i, j);
}
}
return *this;
}
Matrix& operator=(const Type& k){
for(int i = 0; i < Row; i++){
for(int j = 0; j < Col; j++){
(*this)(i, j) = k;
}
}
return *this;
}
まずは代入ですが代入だけで３つもありますね。一番上は一般的な行列同士の代入です。２つ目は分かりにくいのですがA = B + Cのような計算を今後実装するときに必要なモノです。任意の型が入れられるので先ほどのような計算で一時オブジェクト型を入れることができます。最後は行列を任意のスカラーで埋めるための代入です。すべての要素を同じ値にしたいときはこれを使います。
Matrix& operator+(){
return *this;
}
Matrix& operator-(){
for(int i = 0; i < Row; i++){
for(int j = 0; j < Col; j++){
(*this)(i, j) = -(*this)(i, j);
}
}
return *this;
}
お次は符号です。何の変哲もないただの符号です。+3や-5の+とか-のことです。なので+なら行列をそのまま返すし-ならすべての要素を-倍して返します。
Matrix& operator+=(const Matrix& A){
for(int i = 0; i < Row; i++){
for(int j = 0; j < Col; j++){
(*this)(i, j) += A(i, j);
}
}
return *this;
}
Matrix& operator+=(const Type k){
for(int i = 0; i < Row; i++){
for(int j = 0; j < Col; j++){
(*this)(i, j) += k;
}
}
return *this;
}
代入演算子には種類があるのでそれらを実装しておきます。ここでは+=を示していますがほかの演算子もやることは同じです。それぞれの要素同士また要素とスカラーで+=を実行するだけです。
bool operator==(const Matrix& A){
bool result = true;
for(int i = 0; i < Row; i++){
for(int j = 0; j < Col; j++){
result &= ((*this)(i, j) == A(i, j));
}
}
return result;
}
bool operator!=(const Matrix& A){
return !((*this) == A);
}
さらに等号と不等号です。==はそれぞれの要素に対して==を実行して１つでも異なる部分があればfalseが返されるというものです。!=はこれを利用して==の結果を反転させて返します。
//転置
Matrix& T(){
Matrix<Type, Col, Row> A(0);
for(int i = 0; i < Row; i++){
for(int j = 0; j < Col; j++){
A(j, i) = (*this)(i, j);
}
}
return A;
}
転置も難しいことはありません。行と列を順に入れ替えた行列を返すだけです。
次が少し難しいかもしれません。行列式と逆行列ですね。ただしコードで見ると難しいだけで数学でやった計算を機械的に解かせているだけです。私たちは行列を見てこことここを交換すると計算が楽だなあとか考えながら解けますがアルゴリズムではそうもいかないので地道に計算してもらいます。
//行列式
Type det(){
if(Row != Col) return false;
Type result = 1;
if(Row == 2){
result = (*this)(0, 0) * (*this)(1, 1) - (*this)(0, 1) * (*this)(1, 0);
}else if(Row == 3){
result = (*this)(0, 0) * (*this)(1, 1) * (*this)(2, 2)
+ (*this)(1, 0) * (*this)(2, 1) * (*this)(0, 2)
+ (*this)(2, 0) * (*this)(0, 1) * (*this)(1, 2)
- (*this)(2, 0) * (*this)(1, 1) * (*this)(0, 2)
- (*this)(1, 0) * (*this)(0, 1) * (*this)(2, 2)
- (*this)(0, 0) * (*this)(2, 1) * (*this)(1, 2);
}else{
Matrix<double, Row, Col> A(0);
A = *this;
double buf = 0.0;
Type samp = 0.0;
for(int i = 0; i < Row; i++){
//ピボットの選択
int pivot = i;
for(int j = i + 1; j < Row; j++){
if(abs(A(j, i)) > abs(A(pivot, i)))
pivot = j;
}
//ピボットが0なら行列式は0
if(A(pivot, i) == 0)
return 0;
//行の入れ替え
if(i != pivot){
for(int j = 0; j < Row; j++){
samp = A(i, j);
A(i, j) = A(pivot, j);
A(pivot, j) = samp;
}
result *= -1;
}
//上三角化
for(int j = i + 1; j < Row; j++){
buf = A(j, i) / A(i, i);
for(int k = i; k < Row; k++)
A(j, k) -= A(i, k) * buf;
}
result *= A(i, i);
}
}
return result;
}
まず行列式ですがサイズが２の場合は斜め同士をかけて引くおなじみのやり方が簡単です。またサイズ３の場合もさらすの方法で解くことができます。問題は４以上なのですがこれは上三角行列を作ります。そうするとあとは体格成分をかけた値がそのまま行列式の値となります。
//逆行列
Matrix<double, Row, Col> inv(){
if(Row != Col) throw std::runtime_error("matrix size does not match.");
//行列式が0なら逆行列は定義されない
if((*this).det() == 0)
throw std::runtime_error("inverse matrix is not defined.");
Matrix<double, Row, Col> result(0);
Matrix<double, Row, Col * 2> expansion(0);
double samp = 0.0;
double buf = 0.0;
//左にもとの行列を代入
for(int i = 0; i < Row; i++){
for(int j = 0; j < Col; j++){
expansion(i, j) = (*this)(i, j);
}
}
//右に単位行列を代入
for(int i = Col; i < Col * 2; i++){
expansion(i - Row, i) = 1;
}
for(int i = 0; i < Row; i++){
//ピボットの選択
int pivot = i;
for(int j = i + 1; j < Row; j++){
if(abs(expansion(j, i)) > abs(expansion(pivot, i)))
pivot = j;
}
//行の入れ替え
if(i != pivot){
for(int j = 0; j < Col * 2; j++){
samp = expansion(i, j);
expansion(i, j) = expansion(pivot, j);
expansion(pivot, j) = samp;
}
}
//左を単位行列に変換
samp = expansion(i, i);
for(int j = 0; j < Col * 2; j++){
expansion(i, j) /= samp;
}
for(int j = 0; j < Row; j++){
if(j != i){
buf = expansion(j, i);
for(int k = 0; k < Col * 2; k++){
expansion(j, k) -= buf * expansion(i, k);
}
}
}
}
//結果の抽出
for(int i = 0; i < Row; i++){
for(int j = Col; j < Col * 2; j++){
result(i, j - Col) = expansion(i, j);
}
}
return result;
}
逆行列ですがこれはすべてのサイズに対して掃き出し法を適用しています。
size_t row(){
return Row;
}
size_t col(){
return Col;
}
Vector<int, 2> size(){
Vector<int, 2> S({Row, Col});
return S;
}
void print() const {
for (int i = 0; i < Row; ++i) {
matrix[i].print();
}
}
最後に便利そうなのでサイズを返す関数と行列の表示を追加して完了としました。
6.四則演算と式テンプレート
次は四則演算ですね。これには式テンプレートを使うのですがまずはコードをお見せしましょう。
/*～加算～*/
//Matrix + Matrix
template<typename L, typename R, int Row, int Col>
class Add_Mat{
private:
const Matrix<L, Row, Col>& _l;
const Matrix<R, Row, Col>& _r;
public:
Add_Mat(const Matrix<L, Row, Col>& l, const Matrix<R, Row, Col>& r):_l(l), _r(r){}
double operator()(size_t i, size_t j)const{
return _l(i, j) + _r(i, j);
}
};
template<typename L, typename R, int Row, int Col>
inline Add_Mat<L, R, Row, Col> operator+(const Matrix<L, Row, Col>& l, const Matrix<R, Row, Col>& r){
return Add_Mat<L, R, Row, Col>(l, r);
}
//Matrix + Scalar
template<typename L, typename R, int Row, int Col>
class Add_MatSca{
private:
const Matrix<L, Row, Col>& _l;
const R& _r;
public:
Add_MatSca(const Matrix<L, Row, Col>& l, const R& r):_l(l), _r(r){}
double operator()(size_t i, size_t j)const{
return _l(i, j) + _r;
}
};
template<typename L, typename R, int Row, int Col>
inline Add_MatSca<L, R, Row, Col> operator+(const Matrix<L, Row, Col>& l, const R& r){
return Add_MatSca<L, R, Row, Col>(l, r);
}
template<typename L, typename R, int Row, int Col>
inline Add_MatSca<L, R, Row, Col> operator+(const L& l, const Matrix<R, Row, Col>& r){
return Add_MatSca<L, R, Row, Col>(r, l);
}
長いですね。上に書いているのは行列同士の計算です。ただしC = A + Bのような計算をすると一時オブジェクトが発生します。小さな行列ならよいのですがデータのコピーが大きいと無視できないわけですね。そこで式テンプレートを使うことでコピーが必要なくなるという仕組みです。A + BをするときにAdd_Mat(A, B)という関数が呼ばれます。するとここでCに関数が代入されようとします。すると=の定義によって要素ごとに代入しようとします。そこでメンバにあるoperator()が反応して要素ごとの足し算が行われるのでその値がCへと格納されます。これにより一時オブジェクトが発生しません。分かりやすく式で書いておきます。
C = A + B;
C = Add_Mat(A, B);
//代入演算子で各要素ごとに代入しようとする
C(i, j) = Add_Mat(A, B)(i, j);
//メンバ関数で()によって各要素の和が返される
C(i, j) = A(i, j) + B(i, j);
あとはこれをほかの演算子の分も作るだけです。
7.その他の計算
//行列の積
template<typename T, int Row, int Col, int Com>
inline Matrix<T, Row, Col> dot(const Matrix<T, Row, Com>& A, const Matrix<T, Com, Col>& B){
Matrix<T, Row, Col> C(0);
for(int i = 0; i < Row; i++){
for(int j = 0; j < Col; j++){
for(int k = 0; k < Com; k++){
C(i, j) += A(i, k) * B(k, j);
}
}
}
return C;
}
ここまで実装した*=や*は行列の要素同士をかけたのでここでは数学的な行列の積を実装します。
//累乗
template<typename Type, int Size>
Matrix<Type, Size, Size> pow(const Matrix<Type, Size, Size>& A, int n){
Matrix<Type, Size, Size> B(0);
B = A;
for(int i = 1; i < n; i++){
B *= A;
}
return B;
}
最後に何となく実装した累乗の計算です。
8.さいごに
今回は行列の演算を実装しました。とはいっても生成AIにかなり頼った部分はあります。それでも身についたことはたくさんあると信じています。次回がいつになるかは分からないですがニューラルネットワークの実装をしたいですね。その前に多次元の行列を実装するのが先なのでしょうか。まあ気が向いたら書きに来ると思います。
最後に今回のコードのリンクを載せておきます。
Matrix.hpp