• C++ Programming for Financial Engineering
    Highly recommended by thousands of MFE students. Covers essential C++ topics with applications to financial engineering. Learn more Join!
    Python for Finance with Intro to Data Science
    Gain practical understanding of Python to read, understand, and write professional Python code for your first day on the job. Learn more Join!
    An Intuition-Based Options Primer for FE
    Ideal for entry level positions interviews and graduate studies, specializing in options trading arbitrage and options valuation models. Learn more Join!

Designing a matrix header library - using expression templates (ET)

Joined
11/5/14
Messages
294
Points
53
Hi guys,

I would like to design a Matrix class that uses expression templates and avoids temporaries. I would like to start with something very simple, so I decided to just implement matrix addition, first.

I created a class AddExp to represent a compound expression and overloaded the () operator. And the copy-assignment operator Matrix& Matrix::operator=(AddExp& addExp) triggers the evaluation, by the recursively invoking () on each of the elements in this expression.

However, the compiler complains that it can't convert AddExp<Matrix<int,3,3>,Matrix<int,3,3>> to a Matrix<int,3,3>. It puzzles me, why the copy-assignment operator of the Matrix class can't do that.

Note that, if I replace Matrix C by auto C, the code will execute fine, but it's important to be able to write statements of the form Matrix C = <something>.

Compiler Explorer - C++ (x86-64 gcc 12.2)

I would like to ask for some help; if you guys played with expression templates, or find something fundamentally wrong with my code, it would be nice to correct it, and learn from it.

Code:
Compile-time matrix using ET:
#include <array>
#include <iostream>


template <typename LHS, typename RHS, typename T>
class AddExp {
   public:
    AddExp(LHS& lhsExp, RHS& rhsExp) :
          m_lhsExp{lhsExp},
          m_rhsExp{rhsExp} {}

    T operator()(int i, int j) { return m_lhsExp(i, j) + m_rhsExp(i, j); }

    template <typename R>
    AddExp<AddExp<LHS, RHS, T>, R, T> operator+(R& exp) {
        return AddExp<AddExp<LHS, RHS, T>, R, T>(*this, exp);
    }


   private:
    LHS& m_lhsExp;
    RHS& m_rhsExp;
};

template <typename T = double, int ROWS = 3, int COLS = 3>
class Matrix  {
   public:

    Matrix() : m_rows{ROWS}, m_cols{COLS}, data{} {}
    Matrix(std::initializer_list<std::initializer_list<T>> lst) : m_rows{ROWS}, m_cols{COLS}, data{} {
        int i{0};
        for (auto& l : lst) {
            int j{0};
            for (auto& v : l) {
                data[m_rows * i + j] = v;
                ++j;
            }
            ++i;
        }
    }

    T& operator()(int i, int j) { return data[ROWS * i + j]; }

    template <typename RHS>
    AddExp<Matrix<T, ROWS, COLS>, RHS, T> operator+(RHS& exp) {
        return AddExp<Matrix<T, ROWS, COLS>, RHS, T>(*this, exp);
    }

    template <typename RHS>
    Matrix<T, ROWS, COLS>& operator=(const RHS& exp) {
        for (int i{}; i < ROWS; ++i) {
            for (int j{}; j < COLS; ++j) {
                (*this)(i, j) = exp(i, j);
            }
        }

        return (*this);
    }

   private:
    std::array<T, ROWS * COLS> data;
    int m_rows;
    int m_cols;
};

int main() {
    Matrix A{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};

    Matrix B{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};

    Matrix<int, 3, 3> C = A + B;  //Does not compile

    return 0;
}
 
Last edited:
There were two problems above:

  • Initialization v/s Assignment. We need to add a copy-constructor in the Matrix class that accepts a const RHS& exp, because Matrix C = A + B calls the copy constructor when building the object C, and not the assignment operator operator=.
  • const correctness. The overloaded method () should be made const.
The below code now builds and executes as expected.

Lazy Evaluation of coefficient wise matrix operations:
#include <array>
#include <cassert>
#include <iostream>
#include <format>

/**
 * Compile time fixed-size matrix.
 */

template <typename LHS, typename RHS, typename T, int M, int N>
class AddExp;

template <typename LHS, typename RHS, typename T, int M, int N>
class SubExp;

template <typename T, int ROWS, int COLS>
class Matrix;

template <typename T, typename RHS, int ROWS, int COLS>
class ScalarMult
{
    template<typename L, typename R>
    ScalarMult(T k, AddExp<L, R, T, ROWS, COLS>& rhsExp) : scalarConstant{ k }, m_rhsExp{ rhsExp } {}

    template<typename L, typename R>
    ScalarMult(T k, SubExp<L, R, T, ROWS, COLS>& rhsExp) : scalarConstant{ k }, m_rhsExp{ rhsExp } {}

    template<typename R>
    ScalarMult(T k, ScalarMult<T, R, ROWS, COLS>& rhsExp) : scalarConstant{ k }, m_rhsExp{ rhsExp } {}

    ScalarMult(T k, Matrix<T, ROWS, COLS>& rhsExp) : scalarConstant{ k }, m_rhsExp{ rhsExp } {}

    T operator()(int i, int j) const { return scalarConstant * m_rhsExp(i, j); }

private:
    T scalarConstant;
    RHS& m_rhsExp;
};

template <typename LHS, typename RHS, typename T, int M, int N>
class AddExp {
public:
    AddExp(LHS& lhsExp, RHS& rhsExp) : m_lhsExp{ lhsExp }, m_rhsExp{ rhsExp } {}

    T operator()(int i, int j) const { return m_lhsExp(i, j) + m_rhsExp(i, j); }

    template <typename R>
    AddExp<AddExp<LHS, RHS, T, M, N>, R, T, M, N> operator+(R& exp) {
        return AddExp<AddExp<LHS, RHS, T, M, N>, R, T, M, N>(*this, exp);
    }

    template <typename R>
    SubExp<AddExp<LHS, RHS, T, M, N>, R, T, M, N> operator-(R& exp) {
        return SubExp<AddExp<LHS, RHS, T, M, N>, R, T, M, N>(*this, exp);
    }

    template <int COLS>
    Matrix<T, M, COLS> operator*(Matrix<T, N, COLS>& exp) {
        Matrix<T, M, COLS> result;
        if (getCols() != exp.getRows())
            throw std::runtime_error("dimension mismatch!");


        int kMax = exp.getRows();
        for (int i{}; i < M; ++i) {
            for (int j{}; j < COLS; ++j) {
                for (int k{}; k < kMax; ++k) {
                    result(i, j) += (*this)(i, k) * exp(k, j);
                }
                // std::cout << "result(" << i << "," << j << ")" << result(i,
                // j) << "\n";
            }
        }

        return result;
    }

    template <int COLS, typename L, typename R>
    Matrix<T, M, COLS> operator*(const AddExp<L, R, T, N, COLS>& exp2) {
        Matrix<T, M, COLS> result;
        if (getCols() != exp2.getRows())
            throw std::runtime_error("dimension mismatch!");

        int kMax = exp2.getRows();
        for (int i{}; i < M; ++i) {
            for (int j{}; j < COLS; ++j) {
                for (int k{}; k < kMax; ++k) {
                    result(i, j) += (*this)(i, k) * exp2(k, j);
                }
                // std::cout << "result(" << i << "," << j << ")" << result(i,
                // j) << "\n";
            }
        }

        return result;
    }

    LHS& getLHS() { return m_lhsExp; }
    RHS& getRHS() { return m_rhsExp; }
    int getRows() const { return m_lhsExp.getRows(); }
    int getCols() const { return m_lhsExp.getCols(); }

private:
    LHS& m_lhsExp;
    RHS& m_rhsExp;
};

template <typename LHS, typename RHS, typename T, int M, int N>
class SubExp {
public:
    SubExp(LHS& lhsExp, RHS& rhsExp) : m_lhsExp{ lhsExp }, m_rhsExp{ rhsExp } {}

    T operator()(int i, int j) const { return m_lhsExp(i, j) - m_rhsExp(i, j); }

    template <typename R>
    SubExp<SubExp<LHS, RHS, T, M, N>, R, T, M, N> operator-(R& exp) {
        return SubExp<SubExp<LHS, RHS, T, M, N>, RHS, T, M, N>(*this, exp);
    }

    template <typename R>
    AddExp<SubExp<LHS, RHS, T, M, N>, R, T, M, N > operator+(R& exp) {
        return AddExp<SubExp<LHS, RHS, T, M, N>, R, T, M, N>(*this, exp);
    }

    template <int COLS>
    Matrix<T, M, COLS> operator*(Matrix<T, N, COLS>& exp) {
        Matrix<T, M, COLS> result;
        if (getCols() != exp.getRows())
            throw std::runtime_error("dimension mismatch!");


        int kMax = exp.getRows();
        for (int i{}; i < M; ++i) {
            for (int j{}; j < COLS; ++j) {
                for (int k{}; k < kMax; ++k) {
                    result(i, j) += (*this)(i, k) * exp(k, j);
                }
                // std::cout << "result(" << i << "," << j << ")" << result(i,
                // j) << "\n";
            }
        }

        return result;
    }

    template <int COLS, typename L, typename R>
    Matrix<T, M, COLS> operator*(const SubExp<L, R, T, N, COLS>& exp2) {
        Matrix<T, M, COLS> result;
        if (getCols() != exp2.getRows())
            throw std::runtime_error("dimension mismatch!");

        int kMax = exp2.getRows();
        for (int i{}; i < M; ++i) {
            for (int j{}; j < COLS; ++j) {
                for (int k{}; k < kMax; ++k) {
                    result(i, j) += (*this)(i, k) * exp2(k, j);
                }
                // std::cout << "result(" << i << "," << j << ")" << result(i,
                // j) << "\n";
            }
        }

        return result;
    }

    LHS& getLHS() { return m_lhsExp; }
    RHS& getRHS() { return m_rhsExp; }
    int getRows() const { return m_lhsExp.getRows(); }
    int getCols() const { return m_lhsExp.getCols(); }

private:
    LHS& m_lhsExp;
    RHS& m_rhsExp;
};

template <typename T = double, int ROWS = 3, int COLS = 3>
class Matrix {
public:
    Matrix() : m_rows{ ROWS }, m_cols{ COLS }, data{} {}

    Matrix(const Matrix& m)
        : m_rows{ m.m_rows }, m_cols{ m.m_cols }, data{ m.data } {}

    template <typename RHS>
    Matrix(const RHS& exp) : m_rows{ ROWS }, m_cols{ COLS } {
        for (int i{}; i < ROWS; ++i) {
            for (int j{}; j < COLS; ++j) {
                (*this)(i, j) = exp(i, j);
            }
        }
    }

    Matrix(std::initializer_list<std::initializer_list<T>> lst)
        : m_rows{ ROWS }, m_cols{ COLS }, data{} {
        int i{ 0 };
        for (auto& l : lst) {
            int j{ 0 };
            for (auto& v : l) {
                data[m_cols * i + j] = v;
                ++j;
            }
            ++i;
        }
    }

    T& operator()(int i, int j) { return data[COLS * i + j]; }
    T operator()(int i, int j) const { return data[COLS * i + j]; }

    template <typename RHS>
    AddExp<Matrix<T, ROWS, COLS>, RHS, T, ROWS, COLS> operator+(RHS& exp) {
        return AddExp<Matrix<T, ROWS, COLS>, RHS, T, ROWS, COLS>(*this, exp);
    }

    template <typename RHS>
    Matrix<T, ROWS, COLS>& operator=(const RHS& exp) {
        for (int i{}; i < ROWS; ++i) {
            for (int j{}; j < COLS; ++j) {
                (*this)(i, j) = exp(i, j);
            }
        }

        return (*this);
    }

    template <typename L, typename R, int P>
    Matrix<T, ROWS, P> operator*(const AddExp<L, R, T, COLS, P>& exp) {
        Matrix<T, ROWS, P> result;
        if (getCols() != exp.getRows())
            throw std::runtime_error("dimensions mismatch");

        int kMax = exp.getRows();
        for (int i{}; i < ROWS; ++i) {
            for (int j{}; j < P; ++j) {
                for (int k{}; k < kMax; ++k) {
                    result(i, j) += (*this)(i, k) * exp(k, j);
                }
                // std::cout << "result(" << i << "," << j << ")" << result(i,
                // j) << "\n";
            }
        }

        return result;
    }

    template <typename L, typename R, int P>
    Matrix<T, ROWS, COLS> operator*(const SubExp<L, R, T, COLS, P>& exp) {
        Matrix<T, ROWS, P> result;
        if (getCols() != exp.getRows())
            throw std::runtime_error("dimensions mismatch");

        int kMax = exp.getRows();
        for (int i{}; i < m_rows; ++i) {
            for (int j{}; j < P; ++j) {
                for (int k{}; k < kMax; ++k) {
                    result(i, j) += (*this)(i, k) * exp(k, j);
                }
            }
        }

        return result;
    }

    template <int P>
    Matrix<T, ROWS, P> operator*(const Matrix<T, COLS, P>& exp) {
        Matrix<T, ROWS, P> result;
        if (getCols() != exp.getRows())
            throw std::runtime_error("dimensions mismatch");

        int kMax = exp.getRows();
        for (int i{}; i < m_rows; ++i) {
            for (int j{}; j < P; ++j) {
                for (int k{}; k < kMax; ++k) {
                    result(i, j) += (*this)(i, k) * exp(k, j);
                }
                // std::cout << "result(" << i << "," << j << ")" << result(i,
                // j) << "\n";
            }
        }

        return result;
    }

    auto getData() { return data; }
    int getRows() const { return m_rows; }
    int getCols() const { return m_cols; }

private:
    std::array<T, ROWS* COLS> data;
    int m_rows;
    int m_cols;
};


template <typename T, int ROWS, int COLS>
std::ostream& operator<<(std::ostream& stream, Matrix<T, ROWS, COLS>& m) {
    stream << "Matrix<" << m.getRows() << "," << m.getCols() << ">"
        << "[ \n";
    for (int i{ 0 }; i < ROWS; ++i) {
        for (int j{ 0 }; j < COLS; ++j) {
            stream << " " << m(i, j) << " ";
        }
        stream << "\n";
    }
    stream << " ]";

    return stream;
}

template <typename LHS, typename RHS, typename T, int M, int N>
std::ostream& operator<<(std::ostream& stream, const AddExp<LHS, RHS, T, M, N>& m) {
    stream << "Matrix<" << m.getRows() << "," << m.getCols() << ">"
        << "[ \n";
    for (int i{ 0 }; i < m.getRows(); ++i) {
        for (int j{ 0 }; j < m.getCols(); ++j) {
            stream << " " << m(i, j) << " ";
        }
        stream << "\n";
    }
    stream << " ]";

    return stream;
}

using Vector1d = Matrix<double, 1, 1>;
using Vector2d = Matrix<double, 2, 1>;
using Vector3d = Matrix<double, 3, 1>;
using Vector4d = Matrix<double, 4, 1>;


using Matrix2d = Matrix<double, 2, 2>;
using Matrix3d = Matrix<double, 3, 3>;
using Matrix4d = Matrix<double, 4, 4>;

using Vector1f = Matrix<float, 1, 1>;
using Vector2f = Matrix<float, 2, 1>;
using Vector3f = Matrix<float, 3, 1>;
using Vector4f = Matrix<float, 4, 1>;

using Matrix2f = Matrix<float, 2, 2>;
using Matrix3f = Matrix<float, 3, 3>;
using Matrix4f = Matrix<float, 4, 4>;

using Matrix2i = Matrix<int, 2, 2>;
using Matrix3i = Matrix<int, 3, 3>;
using Matrix4i = Matrix<int, 4, 4>;


int main() {
    Matrix<int, 3, 3> A{ {1, 0, 0}, {0, 1, 0}, {0, 0, 1} };

    Matrix<int, 3, 3> B{ {1, 0, 0}, {0, 1, 0}, {0, 0, 1} };

    Matrix<int, 3, 3> C = A + B;

    Matrix<int, 3, 3> D = C * (A + B);

    Matrix<int, 3, 1> x = { {1,2,3} };

    Matrix<int, 3, 1> E = A * x;

    Matrix<int, 3, 3> F = (C + C) * (A + B); 
    std::cout << F;

    return 0;
}
 
Last edited:
I added support for matrix products. Whenever a matrix product is encountered, this is evaluated eagerly, as it is more efficient.

Right now, there's a lot of code duplication - so I hope to refresh my inheritance and polymorphism skills, and create a BinaryExp base class, and have AddExp, SubExp and ScalarMult as its sub-types.
 
Why, oh why?
Hi Daniel, do you mean inheritance with templates complicates design, and makes the code unreadable? Did you mean to suggest, that the design should be kept simple...

I just picked up a couple of C++ books from the library, and would like to apply those ideas, to some small mathematical-finance tasks. e.g. write my own AAD code.
 
No, I mean

1. My feeling is your C++ depth is not sufficient.
2. All this stuff has already been done.
3. It's reinvention of the wheel.
 
No, I mean

1. My feeling is your C++ depth is not sufficient.
2. All this stuff has already been done.
3. It's reinvention of the wheel.
For sure. Out of college, I was only trained on C++ 98. I finished reviewing the first 12 topics from this book - Beginning C++ 20 from novice to professional. I am encountering several new ideas. I am not as proficient, as I'd like to be.

On the bright side, I do refer from time to time, to en.cppreference, or after finishing a small task, ask for code reviews on stackoverflow.

When it comes to the choice of tasks, my time would perhaps be better spent on working on some new problem/taking something forward.

But, working on a new task, and learning language features at the same time sounded to me like a double whammy. So, I'd like to atleast get my programming skills to a green-belt/blue-belt level.
 
I added support for matrix products. Whenever a matrix product is encountered, this is evaluated eagerly, as it is more efficient.

Right now, there's a lot of code duplication - so I hope to refresh my inheritance and polymorphism skills, and create a BinaryExp base class, and have AddExp, SubExp and ScalarMult as its sub-types.
This approach is a bit outdated TBH ... the functional language Haskell can do it in its sleep and has been an influence on C++20.
 
Back
Top