のんびりしているエンジニアの日記

ソフトウェアなどのエンジニア的な何かを書きます。

DenoisingAutoEncoderを実装してみた(C++)

C++ 機械学習
Sponsored Links

皆さんこんにちは
お元気ですか。私は元気です。

さて、今日はAutoEncoderの一種であるDenoisingAutoEncoderを実装してみました。
ublasの機能と補助機能ライブラリを使ってます。作ったのはいいのですが、ぶっちゃけファイル多くて困っています。

理論はこちらを参考にさせて頂きました。

Denoising Autoencodersにおける確率的勾配降下法（数式の導出） - Yusuke Sugomori's Blog

ということでソースコードといきたいところなのですが、記述量が凄まじいので最後に掲載します。

基本的にパラメータをDenoisingAutoEncoderのコンストラクタにセットして、fitを使って学習を行います。

間違いなどがありそうなので、もし見つかったらご連絡してくださいませ。

生身のコードも一応載せておこう…興味が有る人は続きを読むから。

以下のソースコードそんなに分割しなくてもよくない？と思っている人は沢山いらっしゃると思いますが、
今後拡張していく方針なので、あえてこんな形でかいています。もっと拡張するつもりです。

ソースコードをご覧いただけるとわかると思いますが見難いです。多分githubにあげます。

注意事項

Dataset.cpp
適当です。ただただ適当です。

test.cpp
ファイルパスは自分のを使ってください。
今は適当に自分のパス（私の）を設定しています。

疑問点

CrossEntropy誤差関数の部分 1-1　が来たらnanエラーが連発されるのですが、これ解決方法ないでしょうか…。
（今無理＋0.00001とかやってるので…）式はあっているような気はするんだけどなぁ

簡単な解説

DenoisingAutoEncoderの動作の主体は
DenoisingAutoEncoderの持つfitメンバ関数です。

forループ内で行っていることはリコンストラクトまでを実行し、それを元にして
パラメータの更新を行っています。

つまり、リコンストラクト→微分→更新をひたすら繰り返す

はっきりいってコレ以外見る必要がないような…。

void DenoisingAutoEncoder::fit(const ublas::matrix<double> &train_input,const ublas::matrix<double> &train_label,double learning_rate,int training_time){
	ublas::matrix<double> x(train_input);
	for(int i = 0; i < training_time; i++){
		cout << i << "回目" << endl;

		ublas::matrix<double> tilde_x = corruption_input(x, 0.3);
		ublas::matrix<double> y = hidden_value(tilde_x);
		ublas::matrix<double> z = reconstractMatrix(y);

		//リコンストラクト側ののバイアスの微分
		ublas::matrix<double> h_bias_grad = x - z;

		//隠れ層バイアス微分
		ublas::scalar_matrix<double> s_matrix(y.size1(),y.size2(),1.0);
		ublas::matrix<double> right = s_matrix - y;
		ublas::matrix<double> a = prod(h_bias_grad, weight);
		ublas::matrix<double> left = ublas_func::multiMatrixElement(a,y);
		ublas::matrix<double> bias_grad = ublas_func::multiMatrixElement(left, right);

                //weightの微分
		ublas::matrix<double> matrix_grad_1 = prod(trans(tilde_x),bias_grad);
		ublas::matrix<double> matrix_grad_2 = prod(trans(h_bias_grad),y);
		ublas::matrix<double> matrix_grad = matrix_grad_1 + matrix_grad_2;

		//パラメータの更新
		h_bias += learning_rate * ublas_func::mean_by_column(h_bias_grad);
		bias += learning_rate * ublas_func::mean_by_column(bias_grad);
		weight += learning_rate * matrix_grad;

		CrossEntropy *ce = new CrossEntropy();
		cout << "誤差" << ce->execute(x,z) << endl;
	}
}

クラス階層表

f:id:tereka:20140923233057p:plain

簡単な解説（クラスの）

Function:sigmoidやsoftmaxの活性化関数(ActiveFunction)の実装、CrossEntropyなどの実装を行うLossFunction
Layer:隠れ層の実装を想定
Dataset:Mnistのデータ解析
UblasLibrary:便利屋
Util:Ublas以外の便利屋

参考文献

Denoising Autoencodersにおける確率的勾配降下法（数式の導出） - Yusuke Sugomori's Blog

ソースコード

ActiveFunction.hpp

#include <iostream>
#include "Function.hpp"
#include "UblasLibrary.hpp"
#define NDEBUG

enum ActiveFunctionType{
	Sigmoid,
	Softmax
};

class ActiveFunction:public Function{
public:
	virtual ublas::matrix<double> execute(const ublas::matrix<double> &ublas_object) = 0;
};

class SigmoidFunction:public ActiveFunction{
public:
	virtual ublas::matrix<double> execute(const ublas::matrix<double> &ublas_object);
};

class SoftmaxFunction:public ActiveFunction{
public:
	virtual ublas::matrix<double> execute(const ublas::matrix<double> &ublas_object);
};

class ActiveFunctionFuctory{
public:
	ActiveFunction* create(int number);
};

ActiveFunction.cpp

#include <iostream>
#include "ActiveFunction.hpp"

ublas::matrix<double> SigmoidFunction::execute(const ublas::matrix<double> &ublas_object){
	return ublas_func::sigmoid(ublas_object);
}

ublas::matrix<double> SoftmaxFunction::execute(const ublas::matrix<double> &ublas_object){
	return ublas_func::softmax(ublas_object);
}

ActiveFunction* ActiveFunctionFuctory::create(int number){
	if(Sigmoid == number){
		return new SigmoidFunction();
	}
	else{
		return new SoftmaxFunction();
	}
}

Layer.hpp

#pragma once
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/triangular.hpp>
#include <boost/numeric/ublas/lu.hpp> 
#include <time.h>
#include "UblasLibrary.hpp"
#include "ActiveFunction.hpp"
#include "LossFunction.hpp"
#define NDEBUG

using namespace boost::numeric;
using namespace std;

class Layer{
public:
	Layer(int input_n,int input_l);
protected:
	ublas::matrix<double> weight;
	ublas::vector<double> bias; //隠れ層へのバイアス
	ublas::matrix<double> hidden_value(const ublas::matrix<double> &train_input);
};

Layer.cpp

#include "Layer.hpp"

Layer::Layer(const int input_n,const int input_l){
	//思み初期化処理
	//入力ベクトル・層数
	cout << input_n << " " << input_l << endl;
	double a = -4 * sqrt(6. / (input_n + input_l));
	//cout << a << endl;

	weight = ublas_func::generateMatrixRandom(input_n,input_l,-4 * sqrt(6. / (input_n + input_l)),4 * sqrt(6. / (input_n + input_l)));
	bias = ublas_func::generateVectorScalar(input_l, 0.0);
}

ublas::matrix<double> Layer::hidden_value(const ublas::matrix<double> &train_input){
	ActiveFunctionFuctory *acf = new ActiveFunctionFuctory();
	ActiveFunction *af = acf->create(0);
	
	ublas::matrix<double> result = prod(train_input,weight);
	ublas::matrix<double> result2 = ublas_func::sumMatrixAndVector(result,bias);
	return af->execute(result2);
}

AutoEncoder.hpp

#pragma once
#include "Layer.hpp"
#define NDEBUG

class AutoEncoder:public Layer{
public:
	AutoEncoder(int input_n,int input_l);
};

LossFunction.hpp

#include "Function.hpp"
#define NDEBUG

using namespace std;
class LossFunction:public Function{
public:
	virtual double execute(ublas::matrix<double> &t,ublas::matrix<double> &y) = 0;
};

class CrossEntropy:public LossFunction{
public:
	double execute(ublas::matrix<double> &t,ublas::matrix<double> &y);
};

LossFunction.cpp

#include "LossFunction.hpp"

double CrossEntropy::execute(ublas::matrix<double> &t,ublas::matrix<double> &y){
	ublas::scalar_matrix<double> one_matrix(t.size1(),t.size2(),1.0001);
	ublas::scalar_matrix<double> double_one_matrix(t.size1(),t.size2(),0.0001);

	ublas::matrix<double> one_matrix_minus_t = one_matrix - t;
	ublas::matrix<double> one_matrix_minus_y = one_matrix - y;
	ublas::matrix<double> log_one_matrix_minus_y = ublas_func::log(one_matrix_minus_y);

	ublas::matrix<double> up_y = y + double_one_matrix;

	ublas::matrix<double> left_value = ublas_func::multiMatrixElement(t,ublas_func::log(up_y));

	ublas::matrix<double> right_value = ublas_func::multiMatrixElementCheck(one_matrix_minus_t,log_one_matrix_minus_y);	

	ublas::matrix<double> result = left_value + right_value;
	ublas::vector<double> total_vect = ublas_func::sum_by_rows(result);

	double result_number =  - (sum(total_vect) / total_vect.size());
	return result_number;
}

DenoisingAutoEncoder.hpp

#pragma once
#include "AutoEncoder.hpp"
#define NDEBUG

using namespace boost::numeric;

class DenoisingAutoEncoder: public AutoEncoder{
public:
	ublas::vector<double> h_bias; //隠れ層バイアス

	DenoisingAutoEncoder(int input_n,int input_l);
	ublas::matrix<double> corruption_input(const ublas::matrix<double> input,const double corruption_per);
	ublas::matrix<double> reconstractMatrix(const ublas::matrix<double> &input);

	void fit(const ublas::matrix<double> &train_input,double learning_rate,int training_time);
	ublas::matrix<double> reconstract(const ublas::matrix<double> &test_input);
};

DenoisingAutoEncoder.cpp

#include "DenoisingAutoEncoder.hpp"

DenoisingAutoEncoder::DenoisingAutoEncoder(int input_n,int input_l):AutoEncoder(input_n,input_l){
	h_bias = ublas_func::generateVectorScalar(input_n, 0.0);
}

//入力データを欠損させる
ublas::matrix<double> DenoisingAutoEncoder::corruption_input(const ublas::matrix<double> input,const double corruption_per){
	ublas::matrix<double> corrupted_matrix(input);
	for(int i = 0; i < corrupted_matrix.size1(); i++){
		for(int j = 0; j < corrupted_matrix.size2(); j++){
			double d = (double)rand() / RAND_MAX;
			if(d < corruption_per){
				corrupted_matrix(i,j) = 0.0;
			}
		}
	}
	return corrupted_matrix;
}

ublas::matrix<double> DenoisingAutoEncoder::reconstractMatrix(const ublas::matrix<double> &input){
	ActiveFunctionFuctory *acf = new ActiveFunctionFuctory();
	ActiveFunction *af = acf->create(0);
	
	ublas::matrix<double> result = prod(input,trans(weight));
	ublas::matrix<double> result2 = ublas_func::sumMatrixAndVector(result,h_bias);
	return af->execute(result2);
}

void DenoisingAutoEncoder::fit(const ublas::matrix<double> &train_input,double learning_rate,int training_time){
	ublas::matrix<double> x(train_input);
	for(int i = 0; i < training_time; i++){
		clock_t start = clock();

		cout << i << "回目" << endl;

		ublas::matrix<double> tilde_x = corruption_input(x, 0.3);
		ublas::matrix<double> y = hidden_value(tilde_x,0);
		ublas::matrix<double> z = reconstractMatrix(y);

		//りこんすとらくとのバイアスの微分
		ublas::matrix<double> h_bias_grad = x - z;

		//隠れそうバイアス微分
		ublas::scalar_matrix<double> s_matrix(y.size1(),y.size2(),1.0);
		ublas::matrix<double> right = s_matrix - y;
		ublas::matrix<double> a = prod(h_bias_grad, weight);
		ublas::matrix<double> left = ublas_func::multiMatrixElement(a,y);
		ublas::matrix<double> bias_grad = ublas_func::multiMatrixElement(left, right);


		ublas::matrix<double> matrix_grad_1 = prod(trans(tilde_x),bias_grad);
		ublas::matrix<double> matrix_grad_2 = prod(trans(h_bias_grad),y);
		ublas::matrix<double> matrix_grad = matrix_grad_1 + matrix_grad_2;

		//パラメータの更新
		h_bias += learning_rate * ublas_func::mean_by_column(h_bias_grad);
		bias += learning_rate * ublas_func::mean_by_column(bias_grad);
		weight += learning_rate * matrix_grad;

		clock_t loss_calc_start = clock();

		CrossEntropy *ce = new CrossEntropy();
		cout << "誤差" << ce->execute(x,z) << endl;
		//cout << "誤差計算時間 " << (clock() - loss_calc_start) / CLOCKS_PER_SEC << endl;
	}
}

ublas::matrix<double> DenoisingAutoEncoder::reconstract(const ublas::matrix<double> &test_input){
	ublas::matrix<double> y = hidden_value(test_input,0);
	ublas::matrix<double> z = reconstractMatrix(y);

	return z;
}

UblasLibrary.hpp

#pragma once
#include <iostream>
#include <vector>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/triangular.hpp> // 三角行列用のヘッダ．前進消去，後退代入に必要
#include <boost/numeric/ublas/lu.hpp>       // LU分解，前進消去，後退代入用のヘッダ

#include "Util.hpp"
#define NDEBUG

using namespace boost::numeric;

namespace ublas_func{
	/**
	各要素一つ一つに対して、適用する
	関数です。 
	*/
	class Function{
	protected:
		double number; //番号
	public:
		virtual double calc(const double element) = 0;
		void setNumber(const double num){
			number = num;
		}
	};

	class AbsFunction:public Function{
	public:
		virtual double calc(const double element){
			return std::abs(element);
		}
	};

	class LogFunction:public Function{
	public:
		virtual double calc(const double element){
			double result = std::log(element);
			return result;
		}
	};

	class Log10Function:public Function{
	public:
		virtual double calc(const double element){
			return std::log10(element);
		}
	};

	class SqrtFunction:public Function{
	public:
		virtual double calc(const double element){
			return std::sqrt(element);
		}
	};

	class PowerFunction:public Function{
	public:
		virtual double calc(const double element){
			return std::pow(element, number);
		}
	};

	class ExpFunction:public Function{
	public:
		virtual double calc(const double element){
			return std::exp(element);
		}
	};

	class SigmoidFunction:public Function{
	public:
		virtual double calc(const double element){
			return 1.0 / (1.0 + exp(-element));
		}
	};

	class SoftmaxFunction{
	public:
		boost::numeric::ublas::vector<double> calc(const boost::numeric::ublas::vector<double> &ublas_vector){
			double maximum = 0.0;
			for(int i = 0; i < ublas_vector.size(); i++){
				if(ublas_vector[i] > maximum){
					maximum = ublas_vector[i];
				}
			}

			double total = 0.0;
			for(int i = 0; i < ublas_vector.size(); i++){
				total += exp(ublas_vector[i] - maximum);
			}

			boost::numeric::ublas::vector<double> result(ublas_vector.size());
			for(int i = 0; i < ublas_vector.size(); i++){
				result[i] = exp(ublas_vector[i] - maximum) / total;
			}

			return result;
		}
	};

	/** 関数を作る */
	template<typename T>
	BOOST_UBLAS_INLINE T abs(const T &ublas_class){
		Function *fc= new AbsFunction();
		T copy;
		apply_function_to_element(ublas_class,copy,fc);
		return copy;
	}

	template<typename T>
	BOOST_UBLAS_INLINE T log(const T &ublas_class){
		Function *fc = new LogFunction();
		T copy;
		apply_function_to_element(ublas_class,copy,fc);
		return copy;
	}

	template<typename T>
	BOOST_UBLAS_INLINE T log10(const T &ublas_class){
		Function *fc = new Log10Function();
		T copy;
		apply_function_to_element(ublas_class, copy, fc);
		return copy;
	}

	template<typename T>
	BOOST_UBLAS_INLINE T sqrt(const T &ublas_class){
		Function *fc = new SqrtFunction();
		T copy;
		apply_function_to_element(ublas_class,copy, fc);
		return copy;
	}

	template<typename T>
	BOOST_UBLAS_INLINE T pow(const T &ublas_class,double number){
		Function *fc = new PowerFunction();
		fc->setNumber(number);
		T copy;
		apply_function_to_element(ublas_class,copy, fc);
		return copy;
	}

	template<typename T>
	BOOST_UBLAS_INLINE T exp(const T &ublas_class){
		Function *fc = new ExpFunction();
		T copy;
		apply_function_to_element(ublas_class, copy,fc);
		return copy;
	}

	template<typename T>
	BOOST_UBLAS_INLINE T sigmoid(const T &ublas_class){
		Function *fc = new SigmoidFunction();
		T copy;
		apply_function_to_element(ublas_class, copy,fc);
		return copy;
	}

	template<typename T>
	boost::numeric::ublas::matrix<T> softmax(const boost::numeric::ublas::matrix<T> &ublas_class){
		boost::numeric::ublas::matrix<double> result(ublas_class.size1(),ublas_class.size2());
		SoftmaxFunction sf;

		for(int i = 0; i < result.size1(); i++){
			row(result,i) = sf.calc(row(ublas_class,i));
		}
		return result;
	}

	/** vector全てに適用 */
	template<typename T>
	void apply_function_to_element(const boost::numeric::ublas::vector<T> &ublas_vector,boost::numeric::ublas::vector<T> &result_vector,Function *fc){
		int size = ublas_vector.size();
		boost::numeric::ublas::vector<T> result(size);
		for(int i = 0; i < size; i++){
			result[i] = fc->calc(ublas_vector[i]);
		}
	}

	/** matrix全てに適用 */
	template<typename T,typename T2>
	void apply_function_to_element(const boost::numeric::ublas::matrix<T> &ublas_matrix,boost::numeric::ublas::matrix<T2> &result_matrix,Function *fc){
		result_matrix.resize(ublas_matrix.size1(),ublas_matrix.size2());
		for(int i = 0; i < ublas_matrix.size1(); i++){
			for(int j = 0; j < ublas_matrix.size2(); j++){
				result_matrix(i,j) = fc->calc(ublas_matrix(i,j));
			}
		}
	}

	template<typename T,typename T2>
	boost::numeric::ublas::matrix<T> multiMatrixElement(const boost::numeric::ublas::matrix<T> &ublas_matrix1,const boost::numeric::ublas::matrix<T2> &ublas_matrix2){
		boost::numeric::ublas::matrix<T> result(ublas_matrix1.size1(),ublas_matrix1.size2());
		for(int i = 0; i < ublas_matrix1.size1(); i++){
			for(int j = 0; j < ublas_matrix2.size2(); j++){
				result(i,j) = ublas_matrix1(i,j) * ublas_matrix2(i,j);
			}
		}
		return result;
	}

	template<typename T,typename T2>
	boost::numeric::ublas::matrix<T> multiMatrixElementCheck(const boost::numeric::ublas::matrix<T> &ublas_matrix1,const boost::numeric::ublas::matrix<T2> &ublas_matrix2){
		boost::numeric::ublas::matrix<T> result(ublas_matrix1.size1(),ublas_matrix1.size2());
		for(int i = 0; i < ublas_matrix1.size1(); i++){
			for(int j = 0; j < ublas_matrix2.size2(); j++){
				result(i,j) = ublas_matrix1(i,j) * ublas_matrix2(i,j);
				if(std::isnan(result(i,j))){
					cout << result(i,j) << " " << ublas_matrix1(i,j) << " " << ublas_matrix2(i,j) << endl;
				}

				if(std::isinf(result(i,j))){
					cout << result(i,j) << " " << ublas_matrix1(i,j) << " " << ublas_matrix2(i,j) << endl;
				}
			}
		}
		return result;
	}

	/**
	*ここから以下は
	*重要そうな関数を並べています。（逆行列とか）
	*/
	/** 逆行列 */
	template<typename T>
	boost::numeric::ublas::matrix<T> inverse(const boost::numeric::ublas::matrix<T> &ublas_matrix){
		boost::numeric::ublas::matrix<T> copy_matrix(ublas_matrix);

		boost::numeric::ublas::permutation_matrix<> pm(ublas_matrix.size1());
		lu_factorize(copy_matrix,pm);
		boost::numeric::ublas::identity_matrix<double> iden_mat(ublas_matrix.size1());

		lu_substitute(copy_matrix,pm,iden_mat);
		return iden_mat;
	}

	/** 行列式 */
	template<typename T>
	double determinant(const boost::numeric::ublas::matrix<T> &ublas_matrix){
		boost::numeric::ublas::matrix<T> copy_matrix(ublas_matrix);

		boost::numeric::ublas::permutation_matrix<> pm(ublas_matrix.size1());
		lu_factorize(copy_matrix,pm);

		double det = 1.0;

		for(int i = 0; i < pm.size(); i++){
			if(pm(i) != i){
				det *= -1.0;
			}

			det *= copy_matrix(i,i);
		}
		return det;
	}

	/** 擬似逆行列 */
	template<typename T>
	boost::numeric::ublas::matrix<T> psudo(const boost::numeric::ublas::matrix<T> &ublas_matrix){
		return prod(inverse(prod(trace(ublas_matrix),ublas_matrix)),ublas_matrix);
	}

	/**
	ここから以下は、列ごと行ごとに実行するプログラムです。
	Ex:列ごとに平均を取ったベクトルを返す、
	*/
	class FunctionPerVector{
	public:
		virtual double calc(const boost::numeric::ublas::vector<double> &ublas_vector) = 0;
	};

	class MeanFunctionPerVector : public FunctionPerVector{
	public:
		virtual double calc(const boost::numeric::ublas::vector<double> &ublas_vector){
			return sum(ublas_vector) / (double)ublas_vector.size();
		};
	};

	class SumFunctionPerVector : public FunctionPerVector{
	public:
		virtual double calc(const boost::numeric::ublas::vector<double> &ublas_vector){
			return sum(ublas_vector);
		};
	};

	/** 統計 axisは列か行か 0:行 1:列　*/
	template<typename T>
	void apply_to_function_per_vector(const boost::numeric::ublas::matrix<T> &ublas_matrix,boost::numeric::ublas::vector<T> &result_vector,FunctionPerVector *fcpv,int axis){
		if (axis == 0){
			int n_rows = ublas_matrix.size1();
			result_vector.resize(n_rows);
			for(int i = 0; i < n_rows; i++){
				double r = fcpv->calc(row(ublas_matrix,i));
				result_vector[i] = r;
			}
		}else{
			int n_cols = ublas_matrix.size2();
			result_vector.resize(n_cols);
			for(int i = 0; i < n_cols; i++){
				double r = fcpv->calc(column(ublas_matrix,i));
				result_vector[i] = r;
			}
		}
	}

	template<typename T>
	boost::numeric::ublas::vector<T> mean_by_rows(boost::numeric::ublas::matrix<T> &ublas_matrix){
		FunctionPerVector *fcpv = new MeanFunctionPerVector();
		boost::numeric::ublas::vector<T> copy;
		apply_to_function_per_vector(ublas_matrix, copy,fcpv, 0);

		return copy;
	}

	template<typename T>
	boost::numeric::ublas::vector<T> mean_by_column(boost::numeric::ublas::matrix<T> &ublas_matrix){
		FunctionPerVector *fcpv = new MeanFunctionPerVector();
		boost::numeric::ublas::vector<T> copy;
		apply_to_function_per_vector(ublas_matrix, copy,fcpv, 1);

		return copy;
	}

	template<typename T>
	boost::numeric::ublas::vector<T> sum_by_rows(boost::numeric::ublas::matrix<T> &ublas_matrix){
		FunctionPerVector *fcpv = new SumFunctionPerVector();
		boost::numeric::ublas::vector<T> copy;
		apply_to_function_per_vector(ublas_matrix, copy,fcpv, 0);
		return copy;
	}

	template<typename T>
	boost::numeric::ublas::vector<T> sum_by_column(boost::numeric::ublas::matrix<T> &ublas_matrix){
		FunctionPerVector *fcpv = new SumFunctionPerVector();
		boost::numeric::ublas::vector<T> copy;
		apply_to_function_per_vector(ublas_matrix,copy ,fcpv, 1);

		return copy;
	}

	/** こまこまとした計算 */
	template<typename T>
	boost::numeric::ublas::matrix<T> sumMatrixAndVector(const boost::numeric::ublas::matrix<T> &ublas_matrix,const boost::numeric::ublas::vector<T> &ublas_vector){
		boost::numeric::ublas::matrix<T> result(ublas_matrix);

		int n_rows = ublas_matrix.size1();
		for(int i = 0; i < n_rows; i++){
			row(result,i) += ublas_vector; 
		}

		return result;
	}

	/** 以下エラーなどのチェック */
	template<typename T>
	bool checkSquareMatrix(boost::numeric::ublas::matrix<T> &ublas_matrix){
		if(ublas_matrix.size1() == ublas_matrix.size2()){
			return true;
		}else{
			std::cout << "正方行列ではありません。" << std::endl;
			return false;
		}
	}

	template<typename T>
	bool checkVectorSize(const boost::numeric::ublas::vector<T> &ublas_vect1,const boost::numeric::ublas::vector<T> &ublas_vect2){
		if(ublas_vect1.size() == ublas_vect2){
			return true;
		}else{
			std::cout << "ベクトルの長さが一致していません" << std::endl;
			return false;
		}
	}


	/**
	std::vector→ublas::vectorへのconvert 
	*/
	template<typename T>
	boost::numeric::ublas::vector<T> convertToUblasVector(const std::vector<T> &vect1){
		ublas::vector<T> ublas_vect(vect1.size());
		for(int i = 0; i < vect1.size(); i++){
			ublas_vect[i] = vect1[i];
		}
		return ublas_vect;
	}

	template<typename T>
	boost::numeric::ublas::matrix<T> convertToUblasMatrix(const std::vector<std::vector<T> > &vect1){
		boost::numeric::ublas::matrix<double> ublas_matrix(vect1.size(),vect1[0].size());
		for(int i = 0; i < vect1.size(); i++){
			for(int j = 0; j < vect1[i].size(); j++){
				ublas_matrix(i,j) = vect1[i][j];
			}
		}
		return ublas_matrix;
	}

	/** -1 ~ 1までで作成するmatrix　*/
	static boost::numeric::ublas::matrix<double> generateMatrixRandom(int n_row, int n_col,double min = -1.0,double max = 1.0){
		boost::numeric::ublas::matrix<double> result(n_row,n_col);
		for(int i = 0; i < n_row; i++){
			for(int j = 0; j < n_col; j++){
				result(i,j) = rand_area(min, max);
			}
		}
		return result;
	}

	/** -1 ~ 1までで作成するvector　*/
	static boost::numeric::ublas::vector<double> generateVectorRandom(int size,double min = 1,double max = 0){
		boost::numeric::ublas::vector<double> result(size);
		for(int i = 0; i < result.size(); i++){
			result[i] = (double)rand() / (RAND_MAX * 20.0);
		}
		return result;
	}

	/** scalar的なmatrix */
	static boost::numeric::ublas::matrix<double> generateMatrixScalar(int n_row, int n_col,double number){
		boost::numeric::ublas::matrix<double> result(n_row,n_col);
		for(int i = 0; i < n_row; i++){
			for(int j = 0; j < n_col; j++){
				result(i,j) = number;
			}
		}
		return result;
	}

	template<typename T>
	static boost::numeric::ublas::vector<T> generateVectorScalar(int size,T num){
		boost::numeric::ublas::vector<T> result(size);
		for(int i = 0; i < result.size(); i++){
			result[i] = num;
		}
		return result;
	}

	template<typename T>
	int maxIndex(boost::numeric::ublas::vector<T> &ublas_vector){
		int maximumIndex = 0.0;
		double maximum = 0.0;

		for(int i = 0; i < ublas_vector.size(); i++){
			if(ublas_vector[i] > maximum){
				maximumIndex = i;
			}
		}
		return maximumIndex;
	}

	/** 出力 */
	template<typename T>
	void outputSize(boost::numeric::ublas::matrix<T> &ublas_matrix){
		std::cout << ublas_matrix.size1() << " " << ublas_matrix.size2() << std::endl;
	}

	//VectorにMatrixが混じっているかどうか
	template<typename T>
	bool isNanMatrix(boost::numeric::ublas::matrix<T> &ublas_matrix){
		for(int i = 0; i < ublas_matrix.size1(); i++){
			for(int j = 0; j < ublas_matrix.size2(); j++){
				if(std::isnan(ublas_matrix(i,j)) == true){
					return true;
				}
			}
		}
		return false;
	}

	//Matrixにinfが混ざっているか
	template<typename T>
	bool isInfMatrix(boost::numeric::ublas::matrix<T> &ublas_matrix){
		for(int i = 0; i < ublas_matrix.size1(); i++){
			for(int j = 0; j < ublas_matrix.size2(); j++){
				if(std::isinf(ublas_matrix(i,j)) == true){
					return true;
				}
			}
		}
		return false;
	}

	//Vectorにnanが混じっているか
	template<typename T>
	bool isNanVector(boost::numeric::ublas::vector<T> &ublas_vector){
		for(int i = 0; i < ublas_vector.size(); i++){
			if(std::isnan(ublas_vector[i]) == true){
				return true;
			}
		}
		return false;
	}
};

Util.hpp

#include <iostream>
#include <math.h>
#include <random>

using namespace std;

extern double rand_area(double min,double max);

Util.cpp

#include "Util.hpp"

double rand_area(double min,double max){
	random_device rd;
	mt19937 gen(rd());

	std::uniform_real_distribution<double> dis(min, max);
	return dis(gen);
}

test.cpp(main関数)

#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/triangular.hpp>
#include "DenoisingAutoEncoder.hpp"
#include "UblasLibrary.hpp"
#define NDEBUG

using namespace boost::numeric;
using namespace std;

vector<vector<double> > generateTestData(int data_size,int data_length){
	vector<vector<double> > result;

	for(int i = 0; i < data_size; i++){
		vector<double> vect;

		double number = 0.0;
		for(int j = 0; j < data_length; j++){
			if(j % (i + 1) == 0){
				if(number == 0.0){
					number = 1.0;
				}else{
					number = 0.0;
				}
			}
			vect.push_back(number);
		}
		result.push_back(vect);
	}
	return result;
}

void testDenoisingAutoEncoder(){
	int data_length = 60;
	int data_size = 30;

	vector<vector<double> > vect_data = generateTestData(data_size, data_length);
	ublas::matrix<double> input_matrix = ublas_func::convertToUblasMatrix(vect_data);

	DenoisingAutoEncoder da(data_length,30);
	da.fit(input_matrix, label_matrix, 0.1, 100);

	ublas::matrix<double> output = da.reconstract(input_matrix);

	for(int i = 0; i < data_size; i++){
		cout << "入力" << endl;
		cout << row(input_matrix,i) << endl;
		cout << "復元" << endl;
		cout << row(output,i) << endl;
	}
}

int main (void) {
	//testLogistic();
	testDenoisingAutoEncoder();
	//testStackedDenoisingAutoEncoder();
	//testUblas();
}

Makefile

CXX = g++
CFLAGS = -O3 -Wall
CXXFLAGS = -O3 -Wall
SOURCE := $(shell ls *.cpp)
HEADER := $(shell ls *.hpp)
OBJS := $(SOURCE:.cpp=.o)
TARGET = test

all: clean $(TARGET)

$(TARGET):$(OBJS)
	$(CXX) -o $@ $(OBJS)

.c.o:
	$(CC) $(CFLAGS) -c $^
.cpp.o:
	$(CXX) $(CXXFLAGS) -c $^

clean:
	rm -f *.o