[C++]Dlib c++ libraryの基本的な行列計算

# はじめに

　Dlib c++での基本的な操作（主に行列計算）をまとめておきます。私が必要だったと思ったものを書いているので、統一性はありません。Examplesはわかりやすいと思うので、やりたい計算がある場合は、[Examples](http://dlib.net/)で探してみるとよいです。（例えば[行列のexample](http://dlib.net/matrix_ex.cpp.html))。

# 環境

Windows 10 home

clang(7.0.1)

Visual Studio 2017 community

dlib c++ (ver.19.21)

# Dlib c++ libraryとは

　Dlibは、機械学習、最適化計算、画像処理など様々な計算を行うc++ライブラリです。コンパイルに少し時間がかかりますが、ヘッダーオンリーで利用できるため、導入時の苦労が比較的少ないです。利用する場合は、c++11以上のコンパイラを使う必要があります。

# Dlib導入

　[Dlib c++ library](http://dlib.net/)よりdlibライブラリをダウンロードします。ダウンロードファイルを解凍し、解凍フォルダ内のdlibフォルダをコンパイルしたいファイルがあるフォルダーにコピーします。もしくは、環境変数でCPLUS_INCLUDE_PATHを指定しているフォルダにコピーします。

# Dlib での行列計算など

## 行列の宣言と要素の参照

```c++
//#include<dlib/matrix.h>
//行列宣言
dlib::matrix<double, 3, 1> y; //[]
dlib::matrix<double, 3, 3> M;
dlib::matrix<double> M2(3,3);
//初期化
M = 54.2, 7.4, 12.1,
     1, 2, 3,
     5.9, 0.05, 1;
M(0,2) = 10.0;//代入
std::cout << M << std::endl;
/*output
	54.2  7.4  10
	   1    2   3 
	 5.9 0.05   1
*/
```

　基本は、`dlib::matrix<型, 行数, 列数> 変数;`もしくは`dlib::matrix<型> 変数(行数、列数)`で宣言します。要素の参照は、`M(行, 列)`で行います。出力ストリームに`<<`することで、全要素の出力ができます。
## 行列の連結

```c++
//#include <dlib/matrix.h>
dlib::matrix<double, 3, 1> M;
dlib::matrix<double> M_row;
dlib::matrix<double> M_col;
M = 1,2,3;

//行方向の連結
M_row = dlib::join_rows(M, M);
std::cout << M_row;
/*output
1 1
2 2
3 3
*/
//列方向の連結
M_col = dlib::join_cols(M, M);
std::cout << M_col;
/*output
1
2
3
1
2
3
*/
```
行列の連結は、列方向は、`dlib::join_cols`で、行方向は`dlib::join_rows`で行います。

## 行列演算（足し算・引き算・積・逆行列・転置行列）

```c++
//#include <dlib/matrix.h>
dlib::matrix<double, 3, 1> y;
dlib::matrix<double, 3, 1> y1;
dlib::matrix<double, 1, 3> y2;
dlib::matrix<double, 3, 3> M;
dlib::matrix<double, 3, 3> M1;
dlib::matrix<double, 3, 3> M2;
//各要素への掛け算
M2 = M1 * 2;//次元:(3,3)->(3,3)

//足し算
M = M1 + M2;//次元:(3,3)->(3,3)

//引き算
M = M1 - M2;//次元:(3,3)->(3,3)

//行列の積
y = M1 * y1; //次元:(3,3)*(3,1)->(3,1)

//要素の積
M = dlib::pointwise_multiply(M1, M2);//次元(3,3)*.(3,3)->(3,3)

//転置行列
y2 = dlib::trans(y1) * M;//次元(3,1).T*(3,3)->(1,3)

//逆行列
M = M1 * dlib::inv(M1);//次元(3,3)*inv(3,3)->(3,3)(単位行列)

```

　足し算、引き算は`+`、`-`と通常のオペレータを利用します。行列の積は、`*`を利用し、要素の積（アダマール積）は、`dlib::pointwise_multiply`で計算します。転置行列は,`dlib::trans()`、逆行列は`dlib::inv()`を利用します。

## 要素への繰り返し計算

```c++
//#include <dlib/matrix.h>
dlib::matrix<double, 3, 3> M;

for (int i; i < M.nr(); ++i){
    for(int j; j < M.nc(); ++j){
        M(i, j) = 10 * i +j;
    }
}
```

　`for(auto &a:v)`のような書き方はできず、添え字によるループを回す必要があります。行と列の要素数は、`M.nr()`と`M.nc()`で参照することができます

## 連立方程式の解

$$
Ax = b
$$

```c++
//#include <dlib/matrix.h>
//#include <dlib/optimization.h>

dlib::matrix<double, 3, 3> A;
dlib::matrix<double, 3, 1> b;
dlib::matrix<double, 3, 1> x;
dlib::lu_decomposition<dlib::matrix<double>> lu(A);
x = lu.solve(b);
```

　求め方はいろいろあると思いますが、LU分解のクラスである`lu_decomposition`クラスを利用します。コンストラクタで`A`を入力し、`solve`メソッドで`b`を入力します。



## 最小値問題の解

```c++
//#include<dlib/matrix.h>
//#include<dlib/optimization.h>
using vec2 = dlib::matrix<double,2,1>;
//目的関数
auto rosen = [](const vec2 &m){
    double x=m(0), y=m(1);
    return 100.0*pow(y-x*x,2) + pow(1-x,2);
};
auto rosen_derivative = [](const vec2 &m){
    double x=m(0), y=m(1);
    vec2 res;
    res(0) = -400*x*(y - x*x) - 2*(1-x);
    res(1) = 200*(y - x*x);
    return res;
}

//最小値問題
vec2 x = {4, 8};//初期値
dlib::find_min(
    dlib::bfgs_search_strategy(),//最小化の方法
    dlib::objective_delta_stop_strategy(1e-7),//打ち切りの方法
    rosen,//目的関数
    rosen_derivative,//目的関数の勾配
    x,//目的関数の変数
    -1//目的関数値の最大リミット
);
//目的関数の勾配が不明の場合
dlib::find_min_using_approximate_derivatives(
    dlib::bfgs_search_strategy(),//最小値の手法
    dlib::objective_delta_stop_strategy(1e-7) ,//収束のしきい値
    rosen,//目的関数
    x,//目的関数の変数
    -1//目的関数値の最大リミット
);

```

　最小値問題の題材は、DlibのExampleに従い、[Rosenbrock function](https://en.wikipedia.org/wiki/Rosenbrock_function)です。

　目的関数が最小となる引数を`find_min`で求めることができます。目的関数は、上記の例だとラムダ式で与えていますが、関数、関数オブジェクトで与えることも可能です（下方の全コードでは、関数、関数オブジェクトの例も載せています）。`find_min`では目的関数と、目的関数の勾配を与える必要があります。目的関数の戻り値の型が`double`で、目的関数の勾配の戻り値の型が`dlib::matrix`であることが注意点です。
　目的関数の勾配がわからない場合は（もしくは求めるのが面倒な場合は）、`find_min_using_approximate_derivatives`を利用します。
　最小値を求める手法は、引数にオブジェクトで与えます。今回は、`dlib::bfgs_search_strategy`と[BFGS法](https://ja.wikipedia.org/wiki/BFGS%E6%B3%95)を利用しています。他の手法として、[cg法、lbfgs法、newton法](http://dlib.net/optimization.html#find_min)が用意されています。
　収束を決める手法もオブジェクトで与えます。例では、`dlib::objective_delta_stop_strategy`で目的関数の変化が1e-7以下になれば、計算を終了するようにしています。他に勾配ベクトルのノルムで計算を終了する`dlib::gradient_norm_stop_strategy`があります。
　最小値問題と違いはありませんが、最大値問題を解くために、`find_max`も用意されています。



# CSV入出力

```c++
//#include <fstream>
//#include <dlib/matrix.h>

dlib::matrix<double, 3, 3> M;
M = 54.2, 7.4, 12.1,
    1, 2, 3,
    5.9, 0.05, 1;
    
#ファイル出力
std::ofstream ofs("test.csv");
ofs << "a,b,c" << std::endl; //ヘッダー
ofs << dlib::csv << M;
ofs.close();
/* In test.csv
54.2, 7.4, 12.1
1, 2, 3
5.9, 0.05, 1
*/
std::ifstream ifs("test.csv");
std::string tmp;
dilb::matrix<double> M1;
std::getline(ifs, tmp);//1行目のヘッダー読み飛ばし
ifs >> M1;
ifs.close();
```

　csvの書き込みは、`std::ofstream`で行えます。csvらしく、区切りに','をつけたい場合は、`dlib::csv`マニピュレータをつけると実現できます。
　csvの読み込みは、`std::ifstream`でOKです。ヘッダーがある場合は、`getline`などで読み飛ばしましょう。



# 参考
http://dlib.net/

# 以下の実行計算例

```c++:dlib_example.cpp
#include <iostream>
#include <fstream>
#include <dlib/matrix.h>
#include <dlib/optimization.h>
using vec2 = dlib::matrix<double,2,1>;
//最小値問題の目的関数：関数
double rosenFunc(const vec2 &m){
    double x=m(0), y=m(1);
    return 100.0*pow(y-x*x,2) + pow(1-x,2);
}

vec2 rosenDerivativeFunc(const vec2 &m){
    double x=m(0), y=m(1);
    vec2 res;
    res(0) = -400*x*(y - x*x) - 2*(1-x);
    res(1) = 200*(y - x*x);
    return res;
}   
//最小値問題の目的関数：関数オブジェクト
struct RosenFunctor{
    double operator()(const vec2 &m) const {
        double x=m(0), y=m(1);
        return 100.0*pow(y-x*x,2) + pow(1-x,2);
    }
};
struct RosenDerivativeFunctor{
    vec2 operator()(const vec2 &m) const {
        double x=m(0), y=m(1);
        vec2 res;
        res(0) = -400*x*(y - x*x) - 2*(1-x);
        res(1) = 200*(y - x*x);
        return res;
    }
};
int main(){
dlib::matrix<double, 3, 1> y;
dlib::matrix<double, 3, 1> y1;
dlib::matrix<double, 1, 3> y2;
dlib::matrix<double, 3, 3> M;
dlib::matrix<double, 3, 3> M1;
dlib::matrix<double, 3, 3> M2;
dlib::vector<double,3> v;
//初期化
y1 = 3.5,
    1.2,
    7.8;
M1 = 54.2, 7.4, 12.1,
     1, 2, 3,
     5.9, 0.05, 1;
M1(0, 2) = 10.0;

//各要素への掛け算
M2 = M1 * 2;
//足し算
M = M1 + M2;
std::cout << "M1+M2:" << std::endl << M << std::endl;
//引き算
M = M1 - M2;
std::cout << "M1-M2:" << std::endl << M << std::endl;
//行列の積
y = M1 * y1; // 次元(3.3)*(3.1)->(3.1)
std::cout << "M1*y1:" << std::endl << y << std::endl;
//転置行列
y2 = dlib::trans(y1) * M;
std::cout << "y1.T * M:" << std::endl << y2 << std::endl;
//逆行列
M = M1 * dlib::inv(M1);
std::cout << "M1 * M1^-1:" << std::endl << M << std::endl;

//繰り返し計算
std::cout << "loop" << std::endl;
for (int i=0; i < M.nr(); ++i) {
    for(int j=0; j < M.nc(); ++j){
        std::cout << M(i, j) << ", ";
    }
    std::cout << std::endl;
}
std::cout << std::endl;

//連立方程式計算
dlib::matrix<double> A = M1;
dlib::matrix<double> b = y1;
dlib::lu_decomposition<dlib::matrix<double>> lu(A);
std::cout << "Ax=b solution:" << std::endl;
std::cout << lu.solve(b) << std::endl;

//最小値問題
//ラムダ式
auto rosen = [](const vec2 &m){
    double x=m(0), y=m(1);
    return 100.0*pow(y-x*x,2) + pow(1-x,2);
};
auto rosen_derivative = [](const vec2 &m){
    double x=m(0), y=m(1);
    vec2 res;
    res(0) = -400*x*(y - x*x) - 2*(1-x);
    res(1) = 200*(y - x*x);
    return res;
};
vec2 x={4., 8.};
dlib::find_min(
    dlib::bfgs_search_strategy(),
    dlib::objective_delta_stop_strategy(1e-7),
    rosen,
    rosen_derivative,
    x,
    -1
);
std::cout << "Minimization:" << std::endl;
std::cout << "Objective Lambda expression" << std::endl;
std::cout <<"find_min" << std::endl;
std::cout << x << std::endl;

x = {4., 8.};
dlib::find_min_using_approximate_derivatives(
    dlib::bfgs_search_strategy(),//最小値の手法
    dlib::objective_delta_stop_strategy(1e-7) ,//収束のしきい値
    rosen,//最小値を求めたい関数
    x,
    -1);
std::cout << "find_min_using_approximate_derivatives" << std::endl;
std::cout << x << std::endl;

x = {4., 8.};
dlib::find_min(
    dlib::bfgs_search_strategy(),
    dlib::objective_delta_stop_strategy(1e-7),
    rosenFunc,
    rosenDerivativeFunc,
    x,
    -1
);
std::cout << "Function expression:" << std::endl;
std::cout << "find_min"  << std::endl;
std::cout << x << std::endl;

x = {-4., 8.};
dlib::find_min_using_approximate_derivatives(
    dlib::bfgs_search_strategy(),
    dlib::objective_delta_stop_strategy(1e-7),
    rosenFunc,
    x,
    -1
);
std::cout << "find_min_using_approximate_derivatives" << std::endl;
std::cout << x << std::endl;

x = {-4., 8.};
//const RosenFunctor functor;
//const RosenDerivativeFunctor derivative_functor;
dlib::find_min(
    dlib::bfgs_search_strategy(),
    dlib::objective_delta_stop_strategy(1e-7),
    RosenFunctor(),
    RosenDerivativeFunctor(),
    x,
    -1
);
std::cout << "Functor expression:" << std::endl;
std::cout << "find_min" << std::endl;
std::cout << x << std::endl;

x = {-4., 8.} ;
dlib::find_min_using_approximate_derivatives(
    dlib::bfgs_search_strategy(),
    dlib::objective_delta_stop_strategy(1e-7),
    RosenFunctor(),
    x,
    -1
);
std::cout << "find_min_using_approximate_derivatives" << std::endl;
std::cout << x << std::endl;

//csv output and input
std::cout << "csv output:" << std::endl;
std::cout << M1;
std::ofstream ofs("dlib.csv");
ofs <<"a, b, c" << std::endl;
ofs << dlib::csv << M1;
ofs.close();

std::cout << "csv input:" << std::endl;
std::ifstream ifs("dlib.csv");
std::string tmp;
dlib::matrix<double> M_new;
std::getline(ifs, tmp);//skip first rows
ifs >> M_new;
ifs.close();
std::cout << M_new << std::endl;
return 0;
}
```

コンパイル例

```shell
clang++ dlib_example.cpp -o dlib_example
```