Introduction
An important omission in .NET is numerical computation. There are some VB libraries out there for this purpose, but because I had already converted a few Numerical Recipes routines from C to VB4 years ago, and in order to have tight integration with my Color
class, I decided to upgrade these old routines to VB.NET. It was only a small effort to enhance the class with some basic matrix and vector algebra routines, especially with the possibility of overloading in VB.NET.
Background
Some knowledge about linear algebra and numerical computing comes in handy.
Features
Just as my Color
class, this class is not written for speed, instead I hope it is clear and easy to read, and to extend (hint). Note, however that the routines translated from 'Numerical Recipes in C' are NOT very readable and shouldn't be changed (I just know they work, don't ask me how)! There are 3 basic building blocks: a scalar (single), a vector (1 dimensional array of single) and a matrix (2 dimensional array of single). Some of the features of the class are:
- Singular value decomposition based solver, adapted from 'Numerical Recipes in C'. Can solve over-determined systems of linear equations in a least-squares sense.
- Non-linear optimization based solver, adapted from the Nelder-Mead simplex algorithm in 'Numerical recipes in C'. Can find the parameters that minimizes a user-supplied function. See
ImproveRGBTosRGBTransform
in the Color
class for an example of how to use this function.
- Vector norm
- Element-wise addition, subtraction and scaling of vectors and matrices
- Sum and mean of elements of a vector or columns of a matrix
- Maximum and minimum of the elements of a vector or a matrix
- Transpose of a matrix
- Scalar product of a vector-matrix product
- Matrix, vector-matrix and matrix-vector product
- Extraction of sub-matrix
- Getting and settings row and columns of matrix using a vector
- File IO of matrices and vectors
ToString
for matrices and vectors
I realize that this is by no means a complete BLAS class, I just implemented what I needed. I just hope it is helpful to somebody, maybe as a starting class for more advanced functionality (Please let me now, so I can post enhancements back on the website)
Code example
Here is a piece from the Color
class that actually solves 3 sets of linear equations at once (one for each R, G and B color component), resulting in a set of equations that can be used to transform color triplets from one color space to another, e.g. for calibration purposes. Note that although the set of equations is linear with regards to the unknown coefficients, it can be highly non-linear with regards to the input R, G and B values, depending on the value of iTransformType
which determines the type of polynome in R, G and B that will be used, e.g. the third-order polynome R' = a1 * R + a2 * G + a3 * B + a4 * R * G + a5 * R * B + a6 * G * B + a7 * R * G * B
.
Public Shared Sub ComputeColorSpaceTransform(ByVal sC1(,) As Single, _
ByVal sC2(,) As Single, _
ByRef sTransform(,) As Single, _
ByVal iTransformType As ColorSpaceTransform)
Dim iNrColors, i As Integer
iNrColors = sC1.GetUpperBound(0) + 1
Console.WriteLine("Computing " & iTransformType & _
"-term transform using " & iNrColors & " tristimulus values")
Dim sCNL(iTransformType - 1), sA(iNrColors - 1, _
iTransformType - 1), sC(2) As Single
For i = 0 To iNrColors - 1
Algebra.GetMatrixRow(sC1, sC, i)
ColorTripletToNTuple(sC, sCNL)
Algebra.SetMatrixRow(sA, sCNL, i)
Next i
Console.WriteLine("Design matrix is " & vbNewLine _
& Algebra.ToString(sA))
Dim sX(iTransformType - 1, 2) As Single
Algebra.Solve(sA, sX, sC2)
ReDim sTransform(2, iTransformType - 1)
Algebra.Transpose(sX, sTransform)
End Sub
Improvements in current version
I have removed some bugs, and added the non-linear solver. The Nelder-Mead simplex algorithm is very robust and general, albeit slower than some gradient descent type algorithms. In order to be sure to find a global minimum, it can be restarted with a different starting simplex (a simplex is n + 1 dimensional figure, e.g. for a search in 2 dimensional space, a simplex would be a set of 3 2D points, i.e. a triangle, that encloses the potential solution).
To do
The naming of some of the routines may seem bizarre to mathematicians, I need to dig into my university linear algebra courses and find the proper names for some of the operations on vectors and matrices.