Vectors and Matrices

Vectors, matrices, and linear algebra are important tools in nearly all areas of applied mathematics, and Meta.Numerics offers extensive features in this area.

How do I create a vector or matrix?

For small vectors and matrices whose dimensions are known at compile time, here is some example code that does it in a visually suggestive way:

using Meta.Numerics.Matrices;

ColumnVector v = new ColumnVector(2.0, 0.0, -1.0);
ColumnVector w = new ColumnVector(new double[] {1.0, -0.5, 1.5});
             
SquareMatrix A = new SquareMatrix(new double[,] {
    {1, -2, 3},
    {2, -5, 12},
    {0, 2, -10}
});

For more a more complex method of programmatic instantiation, see the next example.

How do I access individual elements?

You can use one-argument indexers to set and get vector components, and two-argument indexers to set and get matrix entries. Here is some code to create a row vector and a random rectangular matrix using indexers:

using System;

RowVector u = new RowVector(4);
for (int i = 0; i < u.Dimension; i++) u[i] = i;

Random rng = new Random(1);
RectangularMatrix B = new RectangularMatrix(4, 3);
for (int r = 0; r < B.RowCount; r++) {
    for (int c = 0; c < B.ColumnCount; c++) {
        B[r, c] = rng.NextDouble();
    }
}

Notice that, as is conventional for .NET languages, indexes are zero-based. Also notice that, as is conventional for matrix notation, the row index comes before the column index.

How do I print a matrix?

The exact output formatting will depend on your application, of course (web page? WPF application? file on disk?), but here is a very simple method that writes to the console:

public static void PrintMatrix(string name, AnyRectangularMatrix M) {
    Console.WriteLine($"{name}=");
    for (int r = 0; r < M.RowCount; r++) {
        for (int c = 0; c < M.ColumnCount; c++) {
            Console.Write("{0,10:g4}", M[r, c]);
        }
        Console.WriteLine();
    }
}

We will use this method extensively in our examples. Notice that our method accepts the abstract AnyRectangularMatrix type, so that it can print any real-valued matrix or vector type. This is a good principal to use in your own code: have your methods accept the most general type they can do to their job.

Can you do vector and matrix arithmetic?

Absolutely! Here are some sums and products involving vectors and matrices:

// Linear combinations of vectors
PrintMatrix("v + 2 w", v + 2.0 * w);

// Matrix-vector multiplication
PrintMatrix("Av", A * v);

// Matrix multiplication
PrintMatrix("B A", B * A);

If you attempt operations on objects with incompatible dimensions (e.g. multiply a 3 X 3 matrix by a 4-dimensional vector), you will get a DimensionMismatchException.

How do I invert a square matrix?

The simplest way is to just call the Inverse method. Here is some code to compute A's inverse and prove that it works:

SquareMatrix AI = A.Inverse();

// This should print a unit matrix
PrintMatrix("A * AI", A * AI);

If you just want to invert an innocuous matrix for a homework assignment, this method will do the job fine. If you want to write a production system that will solve matrix problems in the fastest, most accurate way possible, you should use a decomposition technique (LU Decomposition, QR Decomposition, Singular Value Decomposition), which Meta.Numerics also supports. Keep in mind that not all matrices are invert-able; the Inverse method will throw a DivideByZeroException if given a non-invert-able matrix.

Can I get a transpose?

Yes. Here is some code to print the transposes of v and B:

PrintMatrix("v^T", v.Transpose);
PrintMatrix("B^T", B.Transpose);

As you will be able to see from the printed output, the transpose of a ColumnVector is a RowVector and the transpose of an N X M rectangular matrix is an M X N rectangular matrix.

What's a read-only matrix?

You might expect that matrix transposition is an O(N²) operation in time and memory that copies all the existing matrix entries into a new, independent matrix, but this is not correct. Instead, it is an extremely fast O(1) operation, which is why we expose it as a property rather than a method. We are able to do this because the returned matrix object is merely a new view into the existing storage of the original matrix object. To prevent this from leading to confusion, the returned matrix is read-only; if you attempt to modify any entry, an InvalidOperationException will be raised. Quite often this is no problem at all because you don't need to modify the returned matrix, so you should be quite pleased to have avoided an unnecessary O(N^2) cost.

If you ever do need to modify the returned matrix, just call Copy on the read-only object to obtain an independent copy. Of course, you will then incur the O(N^2) copying cost.

Meta.Numerics uses this technique to avoid (or delay) copying costs in many parts of its matrix API. For example, when you get transformation matrices from a QR decomposition, or eigenvectors from an eigen-decomposition, they are exposed as properties which rapidly return read-only matrices, which can be subsequently copied if necessary.

How do I find a vector norm?

The vector types all expose Norm methods that compute the Euclidean norm:

Console.WriteLine($"|v| = {v.Norm()}");

By the way, the following code will give you the same answer using a dot-product, although it won't be quite as fast:

Console.WriteLine($"sqrt(v^T v) = {Math.Sqrt(v.Transpose * v)}");

How do I find a matrix norm?

Matrix norms are less commonly used than vector norms, but if you need one, we can compute OneNorm, InfinityNorm, FrobeniumNorm, and MaxNorm for any matrix.

What other kinds of matrices does Meta.Numerics support?

There are specialized classes for diagonal matrices, tri-diagonal matrices, and symmetric matrices. These classes include optimizations and functionality specific to their type. For example, only the SymmetricMatrix class has a CholeskyDecomposition method.

What's missing?

We don't yet offer a complex matrix system as rich as our real matrix system. We don't have specialized classes for symmetric tri-diagonal matrices or banded matrices with more than three bands. We don't have specialized types to handle Vandermonde, Toeplitz, circulant, or Hilbert matrices. Please let us know which of these is most important to you and we will get to work on implementing it.

I need a unit matrix.

No problem. Here's how to get one, and some code to show you that it works as expected:

UnitMatrix I = UnitMatrix.OfDimension(3);
PrintMatrix("IA", I * A);

Since a UnitMatrix can only ever represent a unit matrix, it is read-only. If you want to start with a unit matrix and then alter its elements, just call its ToSquareMatrix method to get a modify-able copy.

Can I test matrices for equality?

Yes. Here's some example code that does it for vectors and matrices:

// Should print False
Console.WriteLine(v == w);

// Should print True
Console.WriteLine(I * A == A);

But this is usually a bad idea, for two reasons. First, the entries are floating point numbers, and it's famously almost never a good idea to make decisions based on exact floating point equality. (The finite precision of floating point operations often causes quantities which should mathematically be equal to have very close, but not exactly equal floating point values.) Second, matrix equality testing is an O(N²) operation, which can hide an whole lot of computing time behind an innocuous == sign.