• Lincpp - Linear algebra in C++
  • Summary
  • Introduction - Example of use
  • Description
  • The Vector class
    • Constructors
    • Methods
    • Functions on alg::Vector
    • Operators
  • The Matrix class
    • Shape
    • Constructors
    • Methods
    • Functions on alg::Matrix
    • Operators
  • Convenience functions
  • Internal types: alg::Row and alg::Column

This library implements functionality and a numerical interface for linear algebra objects and constructs, with aim to make numerical computations in C++ as straightforward to work with as in popular environments such as MATLAB and Numpy.

The main goal is ease of use, readability and numerical efficiency. Since these are somewhat contradictory goals, a compromise has been made. In order to preserve a straightforward mathematical syntax, the "return by value/copy" approach has been used. This way, a complex mathematical expression can be written on the right hand side of an attribution by exploiting operator overloads for algebraic classes.

The idea is to allow for rapid prototyping of numerical computations à la Python and, if and when necessary, the developer can easily locate and rewrite the mission-critical bottlenecks by hand without giving up on the very useful numeric data structures.

Below is a short list of examples of what you can do with the functionality defined here:

• Easily declare and print matrix and vector objects.

#include "linear_algebra.hpp"
#include<iostream>
using alg::Matrix;

int main() {
	Matrix A = { {1, 2}, {3, 4} };
	std::cout << A;
}

Output:

{{ 1 2 }
 { 3 4 }}

• Direct multiplication (and other binary operations) and automatic type handling

using alg::Vector;

Vector v{ { 5, 6 } };
Vector mult = A * v;
std::cout << mult;

Output:

{ 17 39 }

• Unary operations and routines such as inverse, determinant, norm etc.

std::cout << A[0].norm() << "\n"; // Note that this is the norm of row 0 of A
std::cout << A.inverse();         // [TODO:] You can also take `A.norm(2)`

Output:

2.23607

{{ -2 1 }
 { 1.5 -0.5 }}

• Operate on rows and columns by reference.

#include <numeric>

Matrix B(3, 3);
std::iota(B.begin(), B.end(), 1);

std::for_each(B.beginCol(), B.endCol(), [](auto column) {
	column += Vector(3, 10);
	column[1] += 100;
});

std::cout << B;

Output:

{{ 11 12 13 }
 { 114 115 116 }
 { 17 18 19 }}

• Scalar operations.

std::cout << 2*A + 1;

Output:

{{ 3 5 }
 { 7 9 }}

For a full example of application, check this submodule on non-linear optimization. It makes use of these definitions on its core routines.

Below is one function defined in that project. It takes a functional object G that represents a vector function, for instance, the gradient of a multivariable function, $G(x) = \nabla f(x)$, and returns another function that computes its Jacobian at any specific point via a central difference formula:

$$ H_{ij} = \frac{[G(\vec x + eI_j) - G(\vec x - eI_j)]_i + [G(\vec x + eI_i) - G(\vec x - eI_i)]_j}{4e}, $$

where $H = (H_{ij})$ is the resulting Jacobian matrix, $I$ is the identity matrix, and $e$ is the desired precision. Here we call the Jacobian $H$ instead of $J$ because in this project the function is tasked to find the Hessian of a function $f(x)$ from its gradient (or its second derivative from its first derivative, in case of a $x$ being a scalar).

#include <functional>
// notice how the expression within the innermost loop automatically evaluates to a double
std::function<Matrix(Vector)> jacobian(std::function<Vector(Vector)> G, int n, double e = 1e-4)
{
    return [G, n, e](Vector x) {
        Matrix H(n, n);
        Matrix M = e * I(n);
        for (size_t i = 0; i < n; i++)
            for (size_t j = 0; j < n; j++)
                H(i, j) = (G(x + M[j]) - G(x - M[j]))[i] / (4 * e) + (G(x + M[i]) - G(x - M[i]))[j] / (4 * e);
        H = 0.5 * (H + t(H)); // this could be more efficiently but less conveniently done by hand
        return H;
    };
}

The module was designed to provide a natural syntax for mathematical operations on vectors and matrices, much like the Numpy package for the Python language or MATLAB. For instance, this means that the dot product between two vectors $\vec v$ and $\vec u$ should be written with simple use of the * operator for ordinary multiplication:

Vector v{ v1, ..., vn };    // Some different ways to initialize a
Vector u = { u1, ..., un }; // Vector with specific values.
double w = v * u;

Of course, as you can see from the above example, the results of operations have reasonable types.

Note

Therefore, the adopted philosophy of operations is that operations that can return a "lesser" type generally do return a "lesser" type.

For instance, as you can see from the examples above, multiplication between matrices and vectors return Vector types (alg::Vector), and the dot product returns a scalar type (of type double).

Nevertheless, there are instances when you don't want the result of a dot product to be reduced to an object of type alg::Vector, or the result might still be a matrix, for example: $(v_1, \dots, v_m)^T\cdot(u_1,\dots,u_m)$. For this reason, operations between two objects of type alg::Matrix still always return a new object of type alg::Matrix. Operations of the type A.col(j) * A.row(i) (where A is necessarily a square alg::Matrix) also yield an alg::Matrix return type.

The header defines its own namespace, alg.

Vector properties and behavior were encapsulated in a class called Vector. Matrix properties and behavior were encapsulated in a class called Matrix.

Note

Methods and operators return new copies of the objects. Few methods act in-place, mostly insertion and access methods. For instance, the reshape method returns a copy, but the setShape method acts in-place.

Vectors do not distinguish between column or row types. When coupled with other objects like matrices and other vectors, the appropriate behavior is inferred. For example, both A*v and v*A, where A is a Matrix and v is a Vector, give the expected results, as v is conceptually transposed accordingly.

alg::Vector defines an object representing a mathematical vector which behaves in C++ similarly to an 1-D NDarray in Python. Its elements can be accessed by (const or non-const) reference via brackets or at method (vector[i] or vector.at(i)).

Mathematical

• + commutative. When used with:
  • scalars
    • Shifts every element by the scalar. Return type: alg::Vector
  • vectors
    • Sums only if operands are of same length. Return type: alg::Vector
  • matrices
    • Sums only if Matrix is row or column of same length. Return type: alg::Vector
• - if
  • unary:
    • scales each element by -1.
  • binary:
    • same behavior as +.
• * commutative with vectors. Non-commutative with matrices. When used with:
  • scalars
    • Scales every element by the scalar. Return type: alg::Vector
  • vectors
    • If operands of same length, does the dot product and returns the result. Return type: double

    It is commutative because it always treats the left-hand side operand as a row-vector and the right-hand side operand as a column vector.

  • matrices
    • Treats the vector as column and multiplies. May throw if shapes misalign. Return type: alg::Vector
• / non-commutative. Can only be used with:
  • scalars
    • same behavior as *.
• += -= *= /= also avaliable.

Access

• [] Returns a double by reference, the entry at the respective index.

Other operators

• = If lenghts are the same, appropriately copies the contents of the right-hand side to the left-hand side. If used with a matrix, the matrix must be either row or column of same length.
• << Sends a representation of the vector object to the designated stream.

alg::Matrix defines an object representing a matrix in row-major order (which behaves like an array of vector objects although it's implemented with contiguous data). Its elements can be accessed via parentheses (matrix(i, j)). Brackets (matrix[i][j]) notation is possible techincally, but the parentheses method is the direct one. matrix[i] will return a Row object that can be further indexed to extract an element at position j.

alg::Shape is a struct keeping a pair of fields m, n of the number of rows and columns, respectively, of a matrix. It's shorthand used for passing this information packaged together. Like any simple object it can be instantiated with brackets (initializer lists), e.g. {2, 3}.

It also features a comparison operator == for testing equality of shapes.

Mathematical

• + commutative. When used with:
  • scalars
    • Shifts every element by the scalar. Return type: alg::Matrix
  • vectors
    • Sums only if matrix is row/column of same length as the vector. Return type: alg::Vector
  • matrices
    • Sums both matrices if they are of same shape. Return type: alg::Matrix
• - if
  • unary:
    • scales each element by -1.
  • binary:
    • same behavior as +.
• * Only commutes with scalars. When used with:
  • scalars
    • Scales every element by the scalar. Return type: alg::Matrix
  • vectors
    • If operands are of appropriate shapes, does product and returns the resulting vector (of same length as the one being operated upon). Return type: alg::Vector
  • matrices
    • Does matrix multiplication and returns a new matrix of appropriate shape. Return type: alg::Matrix
• / non-commutative. Can only be used with:
  • scalars
    • same behavior as *.
• += -= *= /= also avaliable.

Access

• [] Returns a Row object, the row at the respective index.
• () Performant way to get the behavior of [][]. Example: matrix(i, j) == matrix[i][j].

Assignment

• = When used with:
  • matrices
    • If shapes match, appropriately copies the contents of the right-hand side matrix to the left-hand side one.
  • vectors
    • If matrix is a either row or column and of same length, appropriately copies the contents of the right-hand side vector to the left-hand side matrix.
  • row- or column-objects (special vectors)
    • As special kinds of vector, the matrix is further restricted to be of the respective kind in order to proceed with the copy, not of just of either kind.

Miscellaneous

• << Appropriately sends matrix objects to the designated stream.

When the user needs to iterate over the rows or columns of the matrix, or simply needs access to a particular line of the matrix, special objects are created that hold references to all elements situated in that particular line, so that operations and modifications can be done in place with STL functions or user-defined routines. Implementation-wise the objects don't really keep data on each element of the line, they are just convenient devices to iterate over those elements in the matrix.

As such, these objects support iteration and the same numerical operations as alg::Vector. However, as to shape restrictions for those operations, they're as restricted as matrices. That is, while for instance an alg::Vector can be on either side of a dot product with an alg::Matrix, and alg::Column can only ever multiply a matrix from the left if the matrix is row-like. They also don't support insertions or removals.

The most straightforward way to obtain such an object is by calling the methods row and col of alg::Matrix, such as in auto row = matrix.row(3);. Another "legal" manner is to iterate over the matrix with beginRow, endRow, beginCol, and endCol.

Following is a list of methods for these objects:

• begin()
• end()
• getShape() (returns the shape of the original matrix)
• nrows() (for alg::Row)
• norm() and alg::norm(<row/column>)
• nrows() (for alg::Column)
• operator[]
• operator<<
• operator= (with Vector, Matrix and same type)
• size()
• to_string()