sinshu / numflat

A numerical computation library for C#
MIT License
22 stars 3 forks source link
clustering dsp fft linear-algebra math matrix multivariate-analysis numerics vector

NumFlat

NumFlat is a numerical computation library written entirely in C#.

The goal of this project is to provide a lightweight package for handling various mathematical and computational tasks, including linear algebra, multivariable analysis, clustering, and signal processing, using only C#.

Overview

NumFlat provides types named Vec<T> and Mat<T> for representing vectors and matrices. These type names are intentionally chosen to avoid confusion with vector and matrix types (like Vector<T>) from the System.Numerics namespace. Various linear algebra-related operations can be performed on these types through operator overloading and extension methods.

Vec<T> and Mat<T> can hold numerical types that implement the INumberBase<T> interface, which was newly added in .NET 7. The primary supported types are float, double, and Complex. Other types can be used as well, but support beyond simple arithmetic operations is not provided.

Vec<double> vec = [1, 2, 3];

Mat<double> mat =
[
    [1, 2, 3],
    [4, 5, 6],
    [7, 8, 9],
];

Console.WriteLine(mat * vec);

Installation

.NET 8 is required.

The NuGet package is available.

Install-Package NumFlat

Most classes are in the NumFlat namespace.

using NumFlat;

Usage

Create a new vector

A new vector can be created by listing elements inside [].

Code

// Create a new vector.
Vec<double> vector = [1, 2, 3];

// Show the vector.
Console.WriteLine(vector);

Output

Vector 3-Double
1
2
3

Create a new vector from IEnumerable<T>

Vectors can also be created from objects that implement IEnumerable<T>. Since the vector itself is an IEnumerable<T>, it is also possible to call LINQ methods on the vector if needed.

Code

// Some enumerable.
var enumerable = Enumerable.Range(0, 10).Select(i => i / 10.0);

// Create a vector from an enumerable.
var vector = enumerable.ToVector();

// Show the vector.
Console.WriteLine(vector);

Output

Vector 10-Double
  0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9

Indexer access for vectors

Elements in a vector can be accessed or modified through the indexer.

Code

// Create a new vector.
var vector = new Vec<double>(3);

// Element-wise access.
vector[0] = 4;
vector[1] = 5;
vector[2] = 6;

// Show the vector.
Console.WriteLine(vector);

Output

Vector 3-Double
4
5
6

Vector arithmetic

Basic operations on vectors are provided through operator overloading and extension methods.

Code

// Some vectors.
Vec<double> x = [1, 2, 3];
Vec<double> y = [4, 5, 6];

// Addition.
var add = x + y;

// Subtraction.
var sub = x - y;

// Multiplication by a scalar.
var ms = x * 3;

// Division by a scalar.
var ds = x / 3;

// Pointwise multiplication.
var pm = x.PointwiseMul(y);

// Pointwise division.
var pd = x.PointwiseDiv(y);

// Dot product.
var dot = x * y;

// Outer product.
var outer = x.Outer(y);

// L2 norm.
var l2Norm = x.Norm();

// L1 norm.
var l1Norm = x.L1Norm();

// Infinity norm.
var infinityNorm = x.InfinityNorm();

// Normalization.
var normalized = x.Normalize();

Subvector

A subvector can be created from a vector. The subvector acts as a view of the original vector, and changes to the subvector will affect the original vector.

Code

// Some vector.
Vec<double> x = [3, 3, 3, 3, 3];

// Create a subvector of the vector.
var sub = x[1..4];

// Modify the subvector.
sub.Fill(100);

// Show the original vector.
Console.WriteLine(x);

Output

Vector 5-Double
  3
100
100
100
  3

Creating matrix

Matrices can be generated from 2D arrays.

Code

// Create a new matrix.
Mat<double> matrix =
[
    [1, 2, 3],
    [4, 5, 6],
    [7, 8, 9],
];

// Show the matrix.
Console.WriteLine(matrix);

Output

Matrix 3x3-Double
1  2  3
4  5  6
7  8  9

Indexer access for matrices

Elements in a matrix can be accessed or modified through the indexer.

Code

// Creat a new matrix.
var matrix = new Mat<double>(3, 3);

// Element-wise access.
for (var row = 0; row < matrix.RowCount; row++)
{
    for (var col = 0; col < matrix.ColCount; col++)
    {
        matrix[row, col] = 10 * row + col;
    }
}

// Show the matrix.
Console.WriteLine(matrix);

Output

Matrix 3x3-Double
 0   1   2
10  11  12
20  21  22

Matrix arithmetic

Basic operations on matrices are provided through operator overloading and extension methods.

Code

// Some matrices.
Mat<double> x =
[
    [1, 2, 3],
    [0, 1, 2],
    [0, 0, 1],
];

Mat<double> y =
[
    [1, 0, 0],
    [2, 1, 0],
    [3, 2, 1],
];

// Addition.
var add = x + y;

// Subtraction.
var sub = x - y;

// Multiplication.
var mul = x * y;

// Multiplication by a scalar.
var ms = x * 3;

// Division by a scalar.
var ds = x / 3;

// Pointwise multiplication.
var pm = x.PointwiseMul(y);

// Pointwise division.
var pd = x.PointwiseDiv(y);

// Transposition.
var transposed = x.Transpose();

// Trace.
var trace = x.Trace();

// Determinant.
var determinant = x.Determinant();

// Rank.
var rank = x.Rank();

// Inverse.
var inverse = x.Inverse();

// Pseudo-inverse.
var pseudoInverse = x.PseudoInverse();

// Frobenius norm.
var frobeniusNorm = x.FrobeniusNorm();

// L1 norm.
var l1Norm = x.L1Norm();

// L2 norm.
var l2Norm = x.L2Norm();

// Infinity norm.
var infinityNorm = x.InfinityNorm();

Submatrix

A submatrix can be created from a matrix. The submatrix acts as a view of the original matrix, and changes to the submatrix will affect the original matrix.

Code

// Creat a new matrix.
var x = new Mat<double>(5, 5);
x.Fill(3);

// Create a submatrix of the matrix.
var sub = x[1..3, 2..4];

// Modify the subvector.
sub.Fill(100);

// Show the original matrix.
Console.WriteLine(x);

Output

Matrix 5x5-Double
3  3    3    3  3
3  3  100  100  3
3  3  100  100  3
3  3    3    3  3
3  3    3    3  3

Matrix as a set of vectors

Views of rows or columns as vectors can be obtained through the Rows or Cols properties. Similar to a submatrix, changes to the view will affect the original matrix. These properties implement IEnumerable<Vec<T>>, allowing for LINQ methods to be called on collections of vectors.

Code

// Some matrix.
Mat<double> x =
[
    [1, 2, 3],
    [4, 5, 6],
    [7, 8, 9],
];

// Create a view of a row of the matrix.
Vec<double> row = x.Rows[0];

// Create a view of a column of the matrix.
Vec<double> col = x.Cols[1];

// Convert a matrix to a row-major jagged array.
var array = x.Rows.Select(row => row.ToArray()).ToArray();

// Enumerate all the elements in column-major order.
var elements = x.Cols.SelectMany(col => col);

// The mean vector of the row vectors.
var rowMean = x.Rows.Mean();

// The covariance matrix of the column vectors.
var colCov = x.Cols.Covariance();

Matrix decomposition

The LU, QR, Cholesky, SVD, EVD, and GEVD for all three major number types float, double, and Complex are available. Below is an example to decompose a matrix using SVD.

Code

// Some matrix.
Mat<double> x =
[
    [1, 2, 3],
    [4, 5, 6],
    [7, 8, 9],
];

// Do SVD.
var decomposition = x.Svd();

// Decomposed matrices.
var s = decomposition.S.ToDiagonalMatrix();
var u = decomposition.U;
var vt = decomposition.VT;

// Reconstruct the matrix.
var reconstructed = u * s * vt;

// Show the reconstructed matrix.
Console.WriteLine(reconstructed);

Output

Matrix 3x3-Double
1  2  3
4  5  6
7  8  9

Solving linear equations

Linear equations like Ax = b can be solved with the matrix decomposition methods above. Use Solve() method to get the solution vector of a given right-hand side vector b.

Code

// Some coefficient matrix.
Mat<double> a =
[
    [1, 2, 3],
    [1, 4, 9],
    [2, 3, 5],
];

// Do LU decomposition.
var lu = a.Lu();

// The right-hand vector.
Vec<double> b = [1, 2, 3];

// Compute the solution vector.
var x = lu.Solve(b);

// Show the solution.
Console.WriteLine(x);

Output

Vector 3-Double
 2.25
-1.75
 0.75

In-place operations

Most of the operations have an in-place version, which directly modifies the target vector or matrix without creating a new one.

Code

// Some vector.
Vec<double> vector = [1, 2, 3];

// This allocates a new vector.
var normalized = vector.Normalize();

// This modifies the original vector.
vector.NormalizeInplace();

// Some matrix.
Mat<double> matrix =
[
    [1, 2, 3],
    [4, 5, 6],
    [7, 8, 9],
];

// In-place methods can also directly modify a part of vector or matrix.
matrix.Rows[0].NormalizeInplace();

Reducing memory allocation

Most of the operations have a non-allocation version, which requires a destination buffer provided by user. Below is an example of matrix addition without heap allocation.

Code

// Some matrices.
Mat<double> x =
[
    [1, 2],
    [3, 4],
];

Mat<double> y =
[
    [5, 6],
    [7, 8],
];

// This allocates a new matrix.
var answer = x + y;

// The buffer to store the answer.
var destination = new Mat<double>(2, 2);

// This is the same as 'x + y', but does not allocate heap.
Mat.Add(x, y, destination);

Multivariate analyses

The NumFlat.MultivariateAnalyses namespace provides functionality related to multivariate analysis. It currently supports the following methods, with plans to add more methods in the future.

// Read some data. IReadOnlyList<Vec> xs = ReadSomeData();

// Do PCA. var pca = xs.Pca();

foreach (var x in xs) { // Transform the vector based on the PCA. var transformed = pca.Transform(x); }

// Read some label. IReadOnlyList ys = ReadSomeLabel();

// Do LDA. var lda = xs.Lda(ys);

foreach (var x in xs) { // Transform the vector based on the LDA. var transformed = lda.Transform(x); }


### Distributions
The `NumFlat.Distributions` namespace provides functionality related to probability distributions.
It currently supports the multivariate Gaussian distribution and its diagonal covariance matrix variation.
Maximum likelihood estimation from data, probability density function calculation, and random number generation are possible.
#### Code
```cs
using NumFlat.Distributions;

// Read some data.
IEnumerable<Vec<double>> xs = ReadSomeData();

// Compute the maximum likelihood Gaussian distribution.
var gaussian = xs.ToGaussian();

foreach (var x in xs)
{
    // Compute the log PDF.
    var pdf = gaussian.LogPdf(x);
}

Clustering

The NumFlat.Clustering namespace provides functionality related to clustering. It currently supports the following methods, with plans to add more methods in the future.

// Read some data. IReadOnlyList<Vec> xs = ReadSomeData();

// Compute a k-means model with 3 clusters. var kMeans = xs.ToKMeans(3);

// Compute a GMM with 3 clusters. var gmm = xs.ToGmm(3);


### Signal processing
The `NumFlat.SignalProcessing` namespace provides functionality related to signal processing.
It currently supports the following methods, with plans to add more methods in the future.
* Extract frames using the window function
* Overlap addition
* FFT and IFFT
* STFT and ISTFT
* Convolution
* Resampling
#### Code
```cs
using NumFlat.SignalProcessing;

// Some complex vector.
var samples = new Vec<Complex>(256);
samples[0] = 1;

// Do FFT.
var spectrum = samples.Fft();

// Do IFFT.
samples = spectrum.Ifft();

Leveraging OpenBLAS

NumFlat's vector and matrix memory formats are compatible with BLAS and LAPACK, allowing you to leverage faster linear algebra libraries like OpenBLAS. The following code example demonstrates how to perform SVD using a C# binding of OpenBLAS, OpenBlasSharp. While NumFlat provides its own SVD implementation, the SVD implementation in OpenBLAS can be expected to offer performance improvements, especially for large matrices.

Code

Mat<double> a =
[
    [1, 2, 3],
    [4, 5, 6],
    [7, 8, 9],
];

var s = new Vec<double>(Math.Min(a.RowCount, a.ColCount));
var u = new Mat<double>(a.RowCount, a.RowCount);
var vt = new Mat<double>(a.ColCount, a.ColCount);
var work = new Vec<double>(s.Count - 1);

fixed (double* pa = a.Memory.Span)
fixed (double* ps = s.Memory.Span)
fixed (double* pu = u.Memory.Span)
fixed (double* pvt = vt.Memory.Span)
fixed (double* pwork = work.Memory.Span)
{
    Lapack.Dgesvd(
        MatrixLayout.ColMajor,
        'A', 'A',
        a.RowCount, a.ColCount,
        pa, a.Stride, // The content of 'a' will be destroyed.
        ps,
        pu, u.Stride,
        pvt, vt.Stride,
        pwork);
}

Console.WriteLine(u * s.ToDiagonalMatrix() * vt);

Output

Matrix 3x3-Double
1  2  3
4  5  6
7  8  9

Todo

License

NumFlat depends on the following libraries.

NumFlat is available under the MIT license.