C++ Matrix Implementation – std::vector vs std::unique_ptr

cmatrix

As part of a hobby project, I needed a rectangular Matrix object to maintain state of a grid. At first, the implementation seemed trivial and unworthy of further discussion: (I haven't included all the code, only the relevant code)

template<typename T>
class Matrix {
    uint64_t rows, columns;
    std::vector<T> _data;

    uint64_t get_flat_index(uint64_t row, uint64_t column) const {
        return row * columns + column;
    }

public:
    Matrix(uint64_t rows, uint64_t columns) :
    rows(rows), columns(columns), _data(rows * columns, {}) {}

    //Auto-generated by compiler
    //Matrix(Matrix const&) = default;
    //Matrix(Matrix &&) = default;
    //Matrix & operator=(Matrix const&) = default;
    //Matrix & operator=(Matrix &&) = default;
    //~Matrix() = default;

    T & operator()(uint64_t row, uint64_t column) {
        return _data[get_flat_index(row, column)];
    }

    T const& operator()(uint64_t row, uint64_t column) const {
        return _data[get_flat_index(row, column)];
    }

    bool is_valid(uint64_t row, uint64_t column) const {
        return row < rows && column < columns;
    }

    T & at(uint64_t row, uint64_t column) {
        if(!is_valid(row, column)) throw std::runtime_error("row/column out of bounds!");
        return operator()(row, column);
    }

    T const& at(uint64_t row, uint64_t column) const {
        if(!is_valid(row, column)) throw std::runtime_error("row/column out of bounds!");
        return operator()(row, column);
    }

    uint64_t get_rows() const {return rows;}
    uint64_t get_columns() const {return columns;}

    void resize(uint64_t new_rows, uint64_t new_columns) {
        if (new_rows == rows && new_columns == columns) return;

        if (new_columns == columns) {
            _data.resize(new_rows * new_columns, {});
        }
        else {
            std::vector<T> new_data(new_rows * new_columns, {});
            for (uint64_t row = 0; row < std::min(rows, new_rows); row++) {
                auto beginning_of_row = _data.begin() + (row * columns);
                auto ending_of_row = beginning_of_row + std::min(columns, new_columns);
                auto beginning_of_new_row = new_data.begin() + (row * new_columns);
                std::copy(beginning_of_row, ending_of_row, beginning_of_new_row);
            }
            _data = std::move(new_data);
        }

        columns = new_columns;
        rows = new_rows;
    }

    //Other code, not related to this post
};

So it all seems pretty great right? I can write stuff like Matrix<int> m(50,50);, m(5, 10) = 17;, try {m.at(52, 47) = 99;} catch (std::runtime_error const& e) {std::cerr << "Whoops!" << std::endl;}, and it all just works, right?

Well, it turns out there's at least one situation where the code misbehaves in a major way:

Matrix<bool> is_tested(60, 60); 
is_tested(30, 40) = true; //Does not compile! Whoops.

Yeah. Turns out that because std::vector<bool> has been specialized, it messes with the integrity of my code.

My initial solution was to write a specialization for Matrix<bool>.

template<>
class Matrix<bool> {
    uint64_t rows, columns;
    std::unique_ptr<bool[]> _data;

    //Duplicated: 
    uint64_t get_flat_index(uint64_t rows, uint64_t columns) {/*...*/}
public:
    Matrix(uint64_t rows, uint64_t columns) :
    rows(rows), columns(columns), _data(std::make_unique<bool[]>(rows * columns)) {}

    //I don't get this for free anymore!
    Matrix(Matrix const& m) : Matrix(m.rows, m.columns) {
        std::copy(m._data.get(), m._data.get() + rows * columns, _data.get());
    }

    //I have to include this manually now.
    Matrix(Matrix &&) = default;

    //More duplicated code...
    bool & operator()(uint64_t row, uint64_t column) {/*...*/}
    bool const& operator()(uint64_t row, uint64_t column) const {/*...*/}
    bool is_valid(uint64_t row, uint64_t column) const {/*...*/}
    bool & at(uint64_t row, uint64_t column) {/*...*/}
    bool const& at(uint64_t row, uint64_t column) const {/*...*/}
    uint64_t get_rows() const {/*...*/}
    uint64_t get_columns() const {/*...*/}

    void resize(uint64_t new_rows, uint64_t new_columns) {
        if (new_rows == rows && new_columns == columns) return;

        std::unique_ptr<bool[]> new_data{ std::make_unique<bool[]>(new_rows * new_columns) };

        if (new_columns == columns) {
            std::copy(
                begin(),
                begin() + ((new_rows < rows) ? new_rows * new_columns : rows * new_columns),
                new_data.get()
            );
        }
        else {
            for (uint64_t row = 0; row < std::min(rows, new_rows); row++) {
                auto beginning_of_row = _data.get() + (row * columns);
                auto ending_of_row = beginning_of_row + std::min(columns, new_columns);
                auto beginning_of_new_row = new_data.get() + (row * new_columns);
                std::copy(beginning_of_row, ending_of_row, beginning_of_new_row);
            }
        }

        _data = std::move(new_data);
        columns = new_columns;
        rows = new_rows;
    }

    //All the other code needs to be duplicated as well!
};

This is, of course, frustrating, not least of which since every time I spot a mistake in one version of the code, I have to fix it in the other, and same goes if I redesign something.

So my next thought was to ditch std::vector<T> entirely, and just specialize around std::unique_ptr<T[]>. This solves the code duplication problem, but it means I can't take advantage of any optimization potential that std::vector<T> offers over std::unique_ptr<T[]>, like smart use of allocators and other benefits, all to ensure that Matrix<bool> works correctly. I tried a version that partitions out the divergent code into a superclass called _matrix_impl<T> that specializes around bool itself, leaving Matrix<T> to not have to specialize anything itself, but there was still a significant amount of code duplication on things like the variable declarations and the get_flat_index code (not to mention a lot of the code not listed here being duplicated) and it created its own nightmare for code maintainability, vis-a-vis inheritance of template superclasses.

So ultimately, my question is: what is the best solution for this situation? Since my code doesn't have things like insert, emplace, or other similar constructs, does it make sense to just use std::unique_ptr<T[]> for everything, since many of the benefits I'd otherwise have access to are moot anyways? If I use std::vector<T>instead, is there a way to gracefully handle Matrix<bool> without dealing with the headache that is std::vector<bool>? Is there a superior third/fourth option I haven't even considered?

Best Answer

Since my code doesn't have things like insert, emplace, or other similar constructs, does it make sense to just use std::unique_ptr<T[]> for everything, since many of the benefits I'd otherwise have access to are moot anyways?

It's not a bad idea from a performance standpoint; given that you matrix is fixed-size, you can get away with just one pointer instead of three, so your matrix objects are going to be slightly lighter weight than when using std::vector as a data backend; also, given that that pointer has now way to be modified outside the constructor, the compiler may be able to be extra smart and avoid re-reading it from your object when performing manipulations intermixed with extern function calls (it's a common cause of slight slowdown with std::vector).

OTOH, you are not getting the copy/assignment stuff for free, if this is important it's for you to judge.

If I use std::vector<T> instead, is there a way to gracefully handle Matrix<bool> without dealing with the headache that is std::vector<bool>?

A possibility that I actually used is to use the std::vector<T>::reference typedefs for your accessors, thus forwarding whatever proxy object std::vector<bool> likes to use straight to your user. So, something like:

typedef std::vector<T>::reference reference;
typedef std::vector<T>::const_reference const_reference;

reference operator()(uint64_t row, uint64_t column) {
    return _data[get_flat_index(row, column)];
}

const_reference operator()(uint64_t row, uint64_t column) const {
    return _data[get_flat_index(row, column)];
}

// ... same with at & co. ...

Incidentally, if you are to implement your Matrix class using std::vector as a backend, you can avoid storing the number of rows - the std::vector already stores the full size, so the height is just one division away (but as usual, check if the size reduction of the Matrix object is worth the extra cost of the division by profiling the code against common scenarios).

Related Solutions

Strassen’s Algorithm – How the Matrix Multiplication Method Was Developed

Apart from Strassen, nobody is able to tell you how Strassen has got his idea. Howeber¹, I can tell you, how you could have found that formula yourself—provided that you are interested in algebraic geometry and representation theory. This also gives you the tools to show that Strassen's formula is as good as it can, or more precisely, that there is no formula computing the product of two 2×2 matrices that uses fewer than 7 multiplications.

Since you are interested by matrices I assume you know basic linear algebra and will be a bit blurry for the more advanced details.

First let be E the set of all linear maps from a plane to a plane. This is basically the set of all 2×2 matrices, but we forget about a particular coordinate system—because, if there were a better coordinate system than the “default one” we could have interest in using it for matrix multiplication. We also denote by E† the dual space of E and by X = P(E⊗E†⊗E†) the projective space associated to the tensor product E⊗E†⊗E†.

An element of X = P(E⊗E†⊗E†) of the special form [c⊗α⊗β] can be interpreted as an elementary operation on matrices, which, in some appopriate coordinate systems, reads a coefficient of a matrix A and a coefficient of a matrix B and writes the product of these coefficients in some matrix C. A general element of X is a combination of these elementary operations, so the product π of two matrices, understood as a map from P(E)×P(E) to P(E), is a point in X.

The usual matrix product formula and Strassen's formula can be expressed as combinations of these linear operations, so let me denote by W₁ the set of these elementary operations [c⊗α⊗β] and let me describe geometrically their combinations.

Let W₂ be the variety of secants of W₁ in X. It is obtained by taking the (closure of the) union of all lines going through two (generic) points of W₁. We can think of a it as of the set of all combinations of two elemetary operations.

Let W₃ be the variety of secant planes of W₁ in X. It is obtained by taking the (closure of the) union of all planes going through three (generic) points of W₁. We can think of a it as of the set of all combinations of three elemetary operations.

Similarly, we define secant varieties for greater indices. Note that these varieties grow larger and larger, that is W₁⊂W₂⊂W₃⊂⋯ Hence the classical matrix product formula shows that the product of matrices is a point of W₈. Actually

PROPOSITION(Strassen) — The product of matrices π lies in W₇.

As far as I know, Strassen did not put things that way, however this is a geometric point of view on this question. This point of view is very useful, because it also lets you prove that Strassen's formula is the best, that is, that π does not lie in W₆. Geometric methods developped here can also be used for a broader range of problems.

I hope, I caught your curiosity. You can go further by reading this article by Landsberg and Manivel:

http://arxiv.org/abs/math/0601097

¹ I will not fix this typo, because I caught a cold.

C++ – Allow Iteration of Internal Vector Without Leaking Implementation

allow iteration without leaking the internals is exactly what the iterator pattern promises. Of course that is mainly theory so here is a practical example:

class AddressBook
{
  using peoples_t = std::vector<People>;
public:
  using iterator = peoples_t::iterator;
  using const_iterator = peoples_t::const_iterator;

  AddressBook();

  iterator begin() { return people.begin(); }
  iterator end() { return people.end(); }
  const_iterator begin() const { return people.begin(); }
  const_iterator end() const { return people.end(); }
  const_iterator cbegin() const { return people.cbegin(); }
  const_iterator cend() const { return people.cend(); }

private:
  peoples_t people;
};

You provide standard begin and end methods, just like sequences in the STL and implement them simply by forwarding to vector's method. This does leak some implementation detail namely that you're returning a vector iterator but no sane client should ever depend on that so it is imo not a concern. I've shown all overloads here but of course you can start by just providing the const version if clients should not be able to change any People entries. Using the standard naming has benefits: anyone reading the code immediately knows it provides 'standard' iteration and as such works with all common algorithms, range based for loops etc.

Best Answer

Related Solutions

Strassen’s Algorithm – How the Matrix Multiplication Method Was Developed

C++ – Allow Iteration of Internal Vector Without Leaking Implementation

Related Topic