Transposed Convolution in Eigen

main.cpp

#include <iostream>
using std::cout; using std::endl;
#include <Eigen>
using namespace Eigen;

double f(const Matrix3d &a, Matrix3d &y) {
	return 1/2. * (a-y).squaredNorm();
}

Matrix3d conv(Matrix4d &x, Matrix2d &w) {
	Matrix3d ret;
	for (int i = 0; i < 3; ++i)
		for (int j = 0; j < 3; ++j)
			ret(i,j) = w.cwiseProduct(x.block<2,2>(i,j)).sum();
	return ret;
}

int main() {
	Matrix4d x = Matrix4d::Random();
	Matrix2d w = Matrix2d::Random();
	Matrix3d y;
	y << 1, 2, 3,
	     4, 5, 4,
	     3, 2, 1;

	cout << "By difference quotient" << endl;
	Matrix4d dfdx_dq;
	double eps = 0.0000001;
	for (int i = 0; i < 4; ++i) {
		for (int j = 0; j < 4; ++j) {
			x(i,j) += eps;
			double f1 = f(conv(x, w), y);
			x(i,j) -= 2*eps;
			double f2 = f(conv(x, w), y);
			x(i,j) += eps;
			dfdx_dq(i,j) = (f1-f2)/(2*eps);
		}
	}
	cout << dfdx_dq << endl << endl << endl;

	cout << "By back prop" << endl;
	Matrix3d dfda = conv(x, w) - y;
	Matrix4d dfdx_bp; dfdx_bp.setZero();
	for (int i = 0; i < 3; ++i) {
		for (int j = 0; j < 3; ++j) {
			dfdx_bp.block<2,2>(i,j).array() += w.array() * dfda(i,j);
		}
	}
	cout << dfdx_bp << endl << endl << endl;

	cout << "Differences" << endl;
	cout << dfdx_dq - dfdx_bp << endl;
}