Matrix class with space-optimized triangular (lower and upper) variants

Question

I've written a class for square numerical matrices (of doubles, in this case), in C++11. Since I intend to use this class later for implementing an \$LU\$ decomposition, I also subclassed the main Matrix class with LowerTriangular and UpperTriangular matrix classes. Taking advantage of the fact that an \$n\$ by \$n\$ triangular matrix only has \$\frac{n(n + 1)}{2}\$ [potentially nonzero] elements, the latter two classes have been space-optimized accordingly. Alas, that optimization introduced some clumsiness and decision-making regarding the interface, as you will see.

Some quick preliminaries:

#include <valarray>
#include <initializer_list>
#include <iterator>
#include <stdexcept>

using std::size_t;
using std::ptrdiff_t;
typedef size_t index_t;

template<typename T>
using nested_init_list = std::initializer_list<std::initializer_list<T>>;

The Matrix class

Here's the synopsis for the main Matrix class. It uses a single (flat) std::valarray<double> to store its elements. Thus, the subscript operator (which in this case is operator()), needs to perform a small amount of logic to get the corresponding element, by using flattenIndex, which converts an index of the form (index_t i, index_t j) to a single index to subscript the internal flat array. You'll see why I made the subscript operators virtual in the next section. The same applies for choosing to make a custom iterator class instead of just using double * and const double *.

/// Representation of a square matrix
class Matrix {
public:
    /// Create an n*n dense matrix with entries set to 0
    explicit Matrix(size_t n) : Matrix(n, n*n) {}
    /// Construct a matrix from a 2-dimensional init list
    Matrix(const nested_init_list<double> &init);

    // Getters:
    inline size_t size() const { return _n; }
    inline const std::valarray<double> &data() { return _data; }

    // Subscript operators (const and non-const):
    virtual double operator()(index_t i, index_t j) const;
    virtual double &operator()(index_t i, index_t j);

    // Iterator stuff:
    class iterator;
    class const_iterator;
    const_iterator begin() const;
    const_iterator end() const;
    iterator begin();
    iterator end();

protected:
    size_t _n;  /// dimension of matrix
    std::valarray<double> _data;  /// flat array containing the data

    // Constructs an n*n Matrix with custom data size (intended for subclasses)
    Matrix(size_t dimension, size_t dataSize) : _n(dimension), _data(0.0, dataSize) {}

    // Throws exception if (i, j) is out-of-bounds (i.e. i >= _n or j >= _n)
    void checkMatrixBounds(index_t i, index_t j) const;

    // Get the appropriate index for the internal flat array
    virtual inline
    index_t flattenIndex(index_t i, index_t j) const { return _n*i + j; };

private:
    class base_iterator;
};

Implementation

Excluding stuff involving iterators, which I'll put in a separate section.

Matrix::Matrix(const nested_init_list<double> &init) : _n(init.size()), _data(_n*_n) {
    index_t i = 0;
    for (const auto &row : init) {
        if (row.size() != _n)
            throw std::invalid_argument("unevenly sized init list for square matrix");
        for (double num : row) _data[i++] = num;
    }
}

double Matrix::operator()(index_t i, index_t j) const {
    checkMatrixBounds(i, j);
    return _data[flattenIndex(i, j)];
}

double &Matrix::operator()(index_t i, index_t j) {
    checkMatrixBounds(i, j);
    return _data[flattenIndex(i, j)];
}

inline
void Matrix::checkMatrixBounds(index_t i, index_t j) const {
    if (i >= _n or j >= _n) throw std::out_of_range("matrix subscript out of range");
}

---

Triangular subclasses

I made a Triangular abstract subclass, from which both LowerTrinagular and UpperTriangular inherit, to avoid some code duplication. Its main feature, as mentioned earlier, is that it only physically contains \$\frac{n(n + 1)}{2}\$ elements. My goal is to maintain the interface of a general matrix as much as possible, so that algorithms dealing with references to matrices don't need to care, in general, if the matrix is triangular or not. In order to be able to subscript the matrix as usual, though, we need to override the subscript operator.

Triangular

/// Abstract Matrix subclass for triangular matrices
class Triangular : public Matrix {
public:
    double operator()(index_t i, index_t j) const override;
    double &operator()(index_t i, index_t j) override;

protected:
    explicit Triangular(size_t n) : Matrix(n, n*(n + 1)/2) {}
    void initFromList(const nested_init_list<double> &init);

    virtual bool isInTriangle(index_t i, index_t j) const = 0;
    virtual size_t rowSize(index_t i) const = 0;
};

This class declares some (pure) virtual methods, whose definition is dependent on whether the matrix is lower/upper triangular, that are used in the subscript operator. This way the subscript operator only needs to be written once (in the Triangular class) and works for both cases.

The initFromList is a method used in the construction of triangular matrices. Since it needs access to virtual methods, it is not made into a Triangular constructor.

Implementation

I've chosen to implement the modified subscript operators as follows. The const version returns by value; it returns the (i, j)-th element of the matrix as expected: the stored value if (i, j) is within the "triangle", and 0 otherwise. So there's no difference from a dense matrix storing all of the zeros.

double Triangular::operator()(index_t i, index_t j) const {
    checkMatrixBounds(i, j);
    if (not isInTriangle(i, j)) return 0;
    return _data[flattenIndex(i, j)];
}

The non-const version was more difficult to decide on; in fact, I'm not sure whether to keep the following approach. Since the non-const version returns a non-const reference to the element, if the requested (i, j)-th element is outside the triangle, an exception is raised (since modifying an element outside of the "triangle" to a non-zero value would make the matrix non-triangular).

double &Triangular::operator()(index_t i, index_t j) {
    checkMatrixBounds(i, j);
    if (not isInTriangle(i, j))
        throw std::out_of_range("cannot write to null side of triangular matrix");
    return _data[flattenIndex(i, j)];
}

Some other stuff:

void Triangular::initFromList(const nested_init_list<double> &init) {
    index_t flattened = 0;

    index_t i = 0;
    for (const auto &initRow : init) {
        size_t initRowSize = initRow.size();
        if (initRowSize < rowSize(i) or initRowSize > _n)
            throw std::invalid_argument("unevenly sized init list");

        index_t j = 0;
        for (double num : initRow) {
            if (not isInTriangle(i, j)) {
                if (num != 0.0)
                    throw std::invalid_argument("null side of triangular matrix should be zero");
            } else _data[flattened++] = num;
            ++j;
        }

        ++i;
    }
}

LowerTriangular and UpperTriangular

These classes just "concretize" the concepts declared by Triangular; most of the heavy lifting is done by the latter, so there isn't much to comment.

class LowerTriangular : public Triangular {
public:
    explicit LowerTriangular(size_t n) : Triangular(n) {}
    LowerTriangular(const nested_init_list<double> &init) : Triangular(init.size()) {
        initFromList(init);
    }

private:
    inline index_t flattenIndex(index_t i, index_t j) const override { return i*(i + 1)/2 + j; }
    inline bool isInTriangle(index_t i, index_t j) const override { return j <= i; }
    inline size_t rowSize(index_t i) const override { return i; }
};


class UpperTriangular : public Triangular {
public:
    explicit UpperTriangular(size_t n) : Triangular(n) {}
    UpperTriangular(const nested_init_list<double> &init) : Triangular(init.size()) {
        initFromList(init);
    }

private:
    inline index_t flattenIndex(index_t i, index_t j) const override { return i*_n - i*(i + 1)/2 + j; }
    inline bool isInTriangle(index_t i, index_t j) const override { return j >= i; }
    inline size_t rowSize(index_t i) const override { return _n - i; }
};

---

Concerns

Among others:

• The slight inconsistency in the interface for triangular vs. dense matrices (raising an exception when requesting a non-const reference to an element outside the triangle)
• Complexity of the overall design

Answer 1

In the end, since this was only needed (for now) for \$LU\$ decomposition, I ended up ditching the lower/upper subclasses and went for the classical approach instead (as suggested by @harold), which is storing the lower and upper matrices of the decomposition in a single matrix, by taking advantage of the fact that \$L\$ is unit triangular:

$$ \begin{pmatrix} u_{1,1} & u_{1,2} & u_{1, 3} & \cdots & u_{1,n}\\ \ell_{2,1} & u_{2,2} & u_{2, 3} & \cdots & u_{2,n}\\ \ell_{3,1} & \ell_{3,2} & u_{3, 3} & \cdots & u_{3,n}\\ \vdots & \vdots & \ddots & \ddots & \vdots\\ \ell_{n,1} & \ell_{n,2} & \cdots & \ell_{n, n-1} & u_{n, n}\\ \end{pmatrix} $$

where \$\ell_{i,j}\$ and \$u_{i,j}\$ are the elements of \$L\$ and \$U\$ respectively.

So I made an \$LU\$ class that internally stores a matrix as mentioned, with read-only interfaces for the lower/upper matrices:

class LUDecomposition {
public:
    /**
     * Compute the LU decomposition of a matrix.
     * @param mat  matrix to decompose
     * @param tol  numerical tolerance
     * @throws SingularMatrixError if mat is singular
     */
    explicit LUDecomposition(const Matrix &mat, double tol = numcomp::DEFAULT_TOL);
    LUDecomposition() = default;

    class L;
    class U;
    L lower() const;
    U upper() const;

    // [...]

private:
    Matrix _mat;  ///< decomposition matrix (internal data storage)

    // [...]
};

Here are the interfaces for LUDecomposition::L and LUDecomposition::U:

class LUDecomposition::L {
    L(const LUDecomposition &luObj) : _luObj(luObj) {}

    double operator()(index_t i, index_t j) {
        if (j > i) throw std::out_of_range();
        if (i == j) return 1;
        return _luObj._mat(i, j);
    }

private:
    const LUDecomposition &_luObj;
}

class LUDecomposition::U {
    U(const LUDecomposition &luObj) : _luObj(luObj) {}

    double operator()(index_t i, index_t j) {
        if (i > j) throw std::out_of_range();
        return _luObj._mat(i, j);
    }

private:
    const LUDecomposition &_luObj;
}

I'm hesitant on whether these specific interfaces should inherit from Matrix in some way — they do offer (a very reduced subset of) Matrix functionality, but they are fairly different in some regards.

This also has the advantage of not having to store the 1's in the lower matrix, as opposed to the more general approach for triangular matrices. The disadvantages are less conceptual separation between the two and the users having to be aware of the LU decomposition's internal representation.

---

As you can guess, the definitions for LUDecomposition::lower and LUDecomposition::upper are just

LUDecomposition::L LUDecomposition::lower() const { return L(*this); }

LUDecomposition::U LUDecomposition::upper() const { return U(*this); }