Higher Algebra with Operator Overloads (Brian Beckman)

Recently I did a Channel 9 interview with Beth Massi where I walked through a Visual Basic program that used Generics and Operator overloads to perform some higher mathematics. I thought I’d follow up with a post explaining the details of exactly what I did.

Operator overloads with Generics enable some beautiful designs for data types in Higher Algebra, a branch of mathematics, sometimes called Abstract Algebra. Consider fields and vector spaces. I’ll show you operator overloads at THREE levels in a single design. First, background: In this context, a “field” is a collection of unspecified objects, closed under two associative operators, + and *, that obey the distributive law.

Closed means that for any a, b, and c in the field, a+b and a*c are in the field. Associative means

a + (b + c) = (a + b) + c

a * (b * c) = (a * b) * c

Distributive means

a * (b + c) = a * b + a * c

(b + c) * a = b * a + c * a

The field must also have two special members: the additive unit 0 and the multiplicative unit 1, where

a + 0 = 0 + a = a

a * 1 = 1 * a = a

and must have for every a, an additive inverse, -a, such that a + -a = -a + a = 0. There must also be a multiplicative inverse for every element except 0, written 1/a, such that a * 1/a = 1/a * a = 1.

Some authors insist on the commutative laws (a + b = b + a, a * b = b * a), too, but we don’t, here. The most common examples of fields are the Rationals, the Real numbers, the Complex numbers, all of which have commutative addition and multiplication; and the Quaternions, which have non-commutative multiplication.

Don’t confuse the mathematical “field” with a “field” in a record, structure or class.

In each instance of a field, we define + and * to do anything we want so long as they obey the associative and distributive laws. A “vector space over a field” is the set of n-tuples or vectors built up from members of the field and closed under an additional linear combination law, written with no operator symbol, or sometimes with *. If v and w are any two vectors, and f and g are any two members of the underlying field, then f v + g w is a vector, and, furthermore,

f (v + w) = f v + f w

(v + w) g = v g + w g

(f + g) v = f v + g v

w (f + g) = w f + w g

Vector spaces are central in physics and simulation. Imagine 6-vectors of real numbers; such things represent particle states in “phase space” in classical mechanics. Imagine 4-vectors of complex numbers; such things appear in quantum mechanics. Similar structures occur all over Quantum Theory and Gravitation

References:

“Vector Calculus, Linear Algebra, and Differential Forms” [VCLADF], by John H. Hubbard, Barbara Burke Hubbard;
“Lie Groups, Lie Algebras, and Some of Their Applications,” by Robert Gilmore
“The Geometry of Physics,” by Theodore Frankel

Applications for vectors of quaternions are not so easy to come by, but, they are perfectly well defined, and, if we do our software design right, they “just work.” Ditto for vectors over a purely symbolic field or over any other kind. It’s not difficult to support field-like algebraic structures with non-associative multiplication within the same software design. Such things include the Cayley numbers or octonions, but, because of non-associativity, they don’t mesh easily with linear algebra, and that’s where we want to go. Stop with the quaternions.

We extend the design to inner-product spaces, in which every vector has a dual, and to linear algebras: sets of linear transformations of vectors, realized as matrices. With just a little code, we build a general, extensible, optimizable library suitable for physics, engineering, and mathematics in any finite-dimensional vector spaces.

We want three layers:

underlying types
field types, generic over the underlying types
linear-algebra types, generic over the field types

The underlying types comprise built-ins like Double, and custom types like Rational, Complex, Quaternion, and Symbol. Underlying types should implement the basic field operations, but can have many more operations for convenience, like optimized division routines. Operator overloading is a no-brainer for the underlying types.

At the top level of the linear-algebra types, operator overloads are again a no-brainer for vector + vector, matrix * vector, vector * matrix, matrix * matrix, vector * vector, and so on.

At the middle layer, the field types abstract and narrow the operations down to the bare minimum so that the linear-algebra types know what and only what to expect. In this layer, operator overloading is not quite a no-brainer, and there are a couple of different ways to go with a design.

The first paragraph showed the field axioms, and these are what we must design into our field types. One way to do that is to impose the field operations by interface:

Public Interface IField(Of U) 'U is the "UNDERLYING" type
    'Get the underlying value.
    Property UValue() As U
    'Generic users of IField type shouldn't be able to tell
    'whether the Underlying type U has copy-reference
    'semantics or copy-value semantics, so let's insist that
    'providers of IField implement a Dup operation.
    Function Dup() As IField(Of U)
    'This is what fields do (would be nice to have a contract
    'for the distributive law). Must use a design like this
    'since Interfaces cannot host Shared functions in general,
    'and operator overloads in particular.
    Function Add(ByVal that As IField(Of U)) As IField(Of U)
    Function Mul(ByVal that As IField(Of U)) As IField(Of U)
    Function AdditiveInverse(ByVal that As IField(Of U)) _
    As IField(Of U)
    Function MultiplicativeInverse(ByVal that As IField(Of U)) _
    As IField(Of U)
    'Every field must have these. Implement them as Shared and
    'even Const if possible.
    ReadOnly Property Zero() As IField(Of U)
    ReadOnly Property Unity() As IField(Of U)
    'The following is an auxiliary operation for the Normed
    'Division Algebras to support Inner Product, or
    '"Sesquilinear mappings." This computes the complex
    'conjugate, the quaternion conjugate, and so on. If the
    'underlying field does not have a natural dual field, then
    'it's probably self-dual and just implement "Dual" as "Dup."
    Function Dual() As IField(Of U)
End Interface

An upside of this design is that both Structure types and Class types can implement the interface. Thus, the linear-algebra classes are independent of differences of copy semantics. Programmers may use Structure types for run-time speed at both the underlying-type level and at the field level.

Here’s an example of a Structure implementing a field of Doubles

Public Structure SFDouble
    Implements IField(Of Double)
...
    Private mValue As Double
    Public Property UValue() As Double Implements IField(Of Double).UValue
        Get
            Return mValue
        End Get
        Set(ByVal value As Double)
            mValue = value
        End Set
    End Property
    Public Sub New(ByVal d As Double)
        mValue = d
    End Sub
...
    Public Function Add(ByVal that As IField(Of Double)) As IField(Of Double) _
    Implements IField(Of Double).Add
        Return New SFDouble(Me.UValue + that.UValue)
    End Function
    Public Function Dual() As IField(Of Double) Implements IField(Of Double).Dual
        Return New SFDouble(Me.UValue)
    End Function
    Public Function Dup() As IField(Of Double) Implements IField(Of Double).Dup
        Return New SFDouble(Me.UValue)
    End Function
    Public Function AdditiveInverse(ByVal that As IField(Of Double)) As IField(Of Double) _
    Implements IField(Of Double).AdditiveInverse
        'Let the underlying math throw exception here if that.UValue == 0
        Return New FDouble(-Me.UValue)
    End Function
    Public Function MultiplicativeInverse(ByVal that As IField(Of Double)) As IField(Of Double) _
    Implements IField(Of Double).MultiplicativeInverse
        'Let the underlying math throw exception here if that.UValue == 0
        Return New FDouble(1 / Me.UValue)
    End Function
    Public Function Mul(ByVal that As IField(Of Double)) As IField(Of Double) _
    Implements IField(Of Double).Mul
        Return New SFDouble(Me.UValue * that.UValue)
    End Function
    Private Shared sUnity = New SFDouble(1.0)
    Private Shared sZero = New SFDouble(0.0)
    Public ReadOnly Property Unity() As IField(Of Double) Implements IField(Of Double).Unity
        Get
            Return sUnity
        End Get
    End Property
    Public ReadOnly Property Zero() As IField(Of Double) Implements IField(Of Double).Zero
        Get
            Return sZero
        End Get
    End Property
...
End Structure

A Big Downside though, for the purposes of this exercise, is no operator overloads at the field level, because operator overloads
are Shared by definition (i.e., static) and interfaces can’t support static methods (the reason is that static virtual methods don’t make sense in .NET, and all methods in an interface are virtual).

But, we want operator overloads at the field level, and we can get them through an alternative base-class design for fields. We lose Structures at the field level, since particular fields must inherit from a base Class, but the underlying types can still be structures.

Public MustInherit Class AField(Of U) 'U is the "UNDERLYING" type
...
    Private mUValue As U
    Public Property UValue() As U
        Get
            Return mUValue
        End Get
        Set(ByVal value As U)
            mUValue = value
        End Set
    End Property
    'Now, our operator overloads just call virtual MustOverride Functions.
    'Magically, these are NON-COMMUTATIVE in general, just as we want.
    Public Shared Operator +(ByVal e1 As AField(Of U), _
                             ByVal e2 As AField(Of U)) As AField(Of U)
        Return e1.Add(e2)
    End Operator
    Public Shared Operator -(ByVal e1 As AField(Of U), _
                             ByVal e2 As AField(Of U)) As AField(Of U)
        Return e1.Add(e2.AdditiveInverse())
    End Operator
    Public Shared Operator *(ByVal e1 As AField(Of U), _
                             ByVal e2 As AField(Of U)) As AField(Of U)
        Return e1.Mul(e2)
    End Operator
    Public Shared Operator /(ByVal e1 As AField(Of U), _
                             ByVal e2 As AField(Of U)) As AField(Of U)
        Return e1.Mul(e2.MultiplicativeInverse())
    End Operator
    Public Shared Operator -(ByVal that As AField(Of U)) As AField(Of U)
        Return that.AdditiveInverse()
    End Operator
    Public Shared Operator Not(ByVal that As AField(Of U)) As AField(Of U)
        Return that.Dual()
    End Operator
    'This is what fields do (would be nice to have a contract for the
    'distributive law).
    MustOverride Function Add(ByVal that As AField(Of U)) As AField(Of U)
    MustOverride Function Mul(ByVal that As AField(Of U)) As AField(Of U)
    MustOverride Function AdditiveInverse() As AField(Of U)
    MustOverride Function MultiplicativeInverse() As AField(Of U)
    'Every field must have these. Implement them as Shared and even Const if
    'possible.
    MustOverride ReadOnly Property Zero() As AField(Of U)
    MustOverride ReadOnly Property Unity() As AField(Of U)
    'Generic users of IField type shouldn't be able to tell
    'whether the Underlying type U has copy-reference
    'semantics or copy-value semantics, so insist that
    'subclasses of AField implement a Dup operation.
    MustOverride Function Dup() As AField(Of U)
    'The following is an auxiliary operation for the Normed
    'Division Algebras to support Inner Product, or
    '"Sesquilinear mappings." This computes the complex
    'conjugate, the quaternion conjugate, and so on. If the
    'underlying field does not have a natural dual field, then
    'it's probably self-dual and just implement "Dual" as "Dup."
    MustOverride Function Dual() As AField(Of U)
End Class

Now we have operator overloads at the field level, but they’re not virtual (they can’t be). What they do is statically dispatch to
virtual ADD and MUL methods, which are just like the ones we had in the interface design for the field level. Nice, eh? Here are the implementations for the field of Complexes, which are implemented in an underlying Structure type:

Public Class AComplex
    Inherits AField(Of Complex)
...

    Public Overrides Function Add(ByVal that As AField(Of Complex)) As AField(Of Complex)
        Return New AComplex(Me.UValue + that.UValue)
    End Function
    Public Overrides Function AdditiveInverse() As AField(Of Complex)
        Return New AComplex(-Me.UValue.R, -Me.UValue.I)
    End Function
    Public Overrides Function Dual() As AField(Of Complex)
        Return New AComplex(Not Me.UValue)
    End Function
    Public Overrides Function Dup() As AField(Of Complex)
        Return New AComplex(Me)
    End Function
    Public Overrides Function MultiplicativeInverse() As AField(Of Complex)
        Return New AComplex(1.0 / Me.UValue)
    End Function
    Public Overrides Function Mul(ByVal that As AField(Of Complex)) As AField(Of Complex)
        Return New AComplex(Me.UValue * that.UValue)
    End Function
...

    Public Overrides ReadOnly Property Unity() As AField(Of Complex)
        Get
            Return New AComplex(1.0, 0.0)
        End Get
    End Property
    Public Overrides ReadOnly Property Zero() As AField(Of Complex)
        Get
            Return New AComplex(0.0, 0.0)
        End Get
    End Property
    Public Overrides Function ToString() As String
        Return UValue.ToString()
    End Function
End Class

And here is its underlying type:

Imports System.Text
Public Structure Complex
    Private mR As Double
    Public Property R() As Double
        Get
            Return mR
        End Get
        Set(ByVal value As Double)
            mR = value
        End Set
    End Property
    Private mI As Double
    Public Property I() As Double
        Get
            Return mI
        End Get
        Set(ByVal value As Double)
            mI = value
        End Set
    End Property
...
    Public Shared Operator -(ByVal that As Complex) As Complex
        Return New Complex(-that.R, -that.I)
    End Operator
    Public Shared Operator +(ByVal c1 As Complex, ByVal c2 As Complex) As Complex
        Dim result = New Complex()
        result.R = c1.R + c2.R
        result.I = c1.I + c2.I
        Return result
    End Operator
    Public Shared Operator *(ByVal c1 As Complex, ByVal c2 As Complex) As Complex
        Dim result = New Complex()
        result.R = (c1.R * c2.R) - (c1.I * c2.I)
        result.I = (c1.R * c2.I) + (c1.I * c2.R)
        Return result
    End Operator
    Public Shared Operator Not(ByVal c1 As Complex) As Complex
        Return New Complex(c1.R, -c1.I)
    End Operator
...

    Public Shared Operator *(ByVal c As Complex, ByVal scalar As Double) As Complex
        Return New Complex(c.R * scalar, c.I * scalar)
    End Operator
    Public Shared Operator /(ByVal c As Complex, ByVal scalar As Double) As Complex
        Return New Complex(c.R / scalar, c.I / scalar)
    End Operator
    Public Shared Operator /(ByVal c1 As Complex, ByVal c2 As Complex) As Complex
        Return c1 * Not c2 / c2.MagnitudeSquared()
    End Operator
...

End Structure

Now, we have an example of a Field generic on an underlying type, How about the Vector space? I wrote one vector space generic
on the interface-style Field design, and another vector space generic on the base-class Field design. Here is the latter one, it’s prettier:

Imports System.Text
'The only reason I must mention the underlying type, U, here, is to pass it to the
'Generic AField interface. I don't use U anywhere explicitly inside this. If I had
'"Monads" in the Generic type language, then I could thread U through this via a
'"bind" operator, very pretty, but dreamland.
Public Class ANVector(Of U, T As {New, AField(Of U)})
    Private mDimension As Integer
    Public ReadOnly Property Dimension() As Integer
        Get
            Return mDimension
        End Get
    End Property
    Private mComponents() As T
    Default Property Component(ByVal i As Integer) As T
        Get
            Component = mComponents(i)
        End Get
        Set(ByVal Value As T)
            mComponents(i) = Value
        End Set
    End Property
...

    Public Sub New(ByVal that As ANVector(Of U, T))
        mDimension = that.Dimension
        ReDim mComponents(Dimension)
        For i = 1 To Dimension
            'Cannot call New with parameters on Generic types ... Sorry)
            ' --> Me.mComponents(i) = New T(that(i))
            Me.mComponents(i) = that(i).Dup()
        Next
    End Sub
    'Add two vectors. PROPERLY, vectors should themselves implement IModule, etc. Some day...
    Public Shared Operator +(ByVal V1 As ANVector(Of U, T), ByVal V2 As ANVector(Of U, T)) _
    As ANVector(Of U, T)
        If V1.Dimension <> V2.Dimension Then
            Throw New VectorSpaceException("Cannot add vectors of different dimensions")
        Else
            Dim result = New ANVector(Of U, T)(V1) 'clone
            For i = 1 To result.Dimension
                result(i) = result(i) + V2(i)''' OPERATOR OVERLOADS BEING USED
            Next
            Return result
        End If
    End Operator
    'Scale a vector
    Public Shared Operator *(ByVal that As ANVector(Of U, T), ByVal scalar As T) _
    As ANVector(Of U, T)
        Dim result = New ANVector(Of U, T)(that)
        For i = 1 To result.Dimension
            result(i) = result(i) * scalar''' OPERATOR OVERLOADS BEING USED
        Next
        Return result
    End Operator
    Public Shared Operator *(ByVal scalar As T, ByVal that As ANVector(Of U, T)) _
    As ANVector(Of U, T)
        Return that * scalar''' OPERATOR OVERLOADS BEING USED
    End Operator
    Public Function Dual()
        Dim result = New ANVector(Of U, T)(Me.Dimension)
        For i = 1 To result.Dimension
            result(i) = Not result(i)''' OPERATOR OVERLOADS BEING USED
        Next
        Return result
    End Function
    Public Shared Operator Not(ByVal V As ANVector(Of U, T)) As ANVector(Of U, T)
        Return V.Dual()
    End Operator
    'Inner product
    Public Shared Operator *(ByVal Left As ANVector(Of U, T), ByVal Right As ANVector(Of U, T)) _
    As T
        If Left.Dimension <> Right.Dimension Then
            Throw New VectorSpaceException _
            ("Cannot compute inner product of vectors with different dimensions")
        Else
            Dim result As New T()
            Dim temp As New T()
            For i = 1 To Left.Dimension
                result = result + Left(i) * Not Right(i)''' OPERATOR OVERLOADS BEING USED
            Next
            Return result
        End If
    End Operator
...
End Class
'================================================================

Imports System.Text
Public Class AMNMatrix(Of U, T As {New, AField(Of U)})
    Private mRowCount As Integer
    Public ReadOnly Property Rows() As Integer
        Get
            Return mRowCount
        End Get
    End Property
    Private mColumnCount As Integer
    Public ReadOnly Property Columns() As Integer
        Get
            Return mColumnCount
        End Get
    End Property
    Private mElements As T(,)
    Default Property Element(ByVal i As Integer, ByVal j As Integer) As T
        Get
            Return mElements(i, j)
        End Get
        Set(ByVal Value As T)
            mElements(i, j) = Value
        End Set
    End Property
...
    'Scale a matrix
    Public Shared Operator *(ByVal scalar As T, ByVal that As AMNMatrix(Of U, T)) _
    As AMNMatrix(Of U, T)
        Dim result = New AMNMatrix(Of U, T)(that)
        For i = 1 To that.Rows
            For j = 1 To that.Columns
                result(i, j) = result(i, j) * scalar ''' OPERATOR OVERLOADS BEING USED
            Next
        Next
        Return result
    End Operator
    'Scale a matrix
    Public Shared Operator *(ByVal that As AMNMatrix(Of U, T), ByVal scalar As T) _
    As AMNMatrix(Of U, T)
        Return scalar * that''' OPERATOR OVERLOADS BEING USED
    End Operator

    'Matrix plus matrix 
    Public Shared Operator +(ByVal M1 As AMNMatrix(Of U, T), ByVal M2 As AMNMatrix(Of U, T)) _
    As AMNMatrix(Of U, T)
        If M1.Rows <> M2.Rows Or M1.Columns <> M2.Columns Then
            Throw New VectorSpaceException("Cannot add matrices of different dimensions")
        Else
            Dim result = New AMNMatrix(Of U, T)(M1)
            For i = 1 To M1.Rows
                For j = 1 To M1.Columns
                    result(i, j) = result(i, j) + M2(i, j) ''' OPERATOR OVERLOADS BEING USED
                Next
            Next
            Return result
        End If
    End Operator
    'Matrix times a vector
    Public Shared Operator *(ByVal m As AMNMatrix(Of U, T), ByVal v As ANVector(Of U, T)) _
    As ANVector(Of U, T)
        If m.Columns <> v.Dimension Then
            Throw New VectorSpaceException("The number of columns in the matrix must equal " & _
                                           "the dimension of the vector")
        Else
            Dim result = New ANVector(Of U, T)(m.Rows)
            For i = 1 To m.Rows
                For j = 1 To m.Columns
                    result(i) = result(i) + v(j) * m(i, j) ''' OPERATOR OVERLOADS BEING USED
                Next
            Next
            Return result
        End If
    End Operator
    'Vector times a Matrix
    Public Shared Operator *(ByVal v As ANVector(Of U, T), ByVal m As AMNMatrix(Of U, T)) _
    As ANVector(Of U, T)
        If m.Rows <> v.Dimension Then
            Throw New VectorSpaceException("The dimension of the vector must equal" & _
                                           "the number of rows in the matrix")
        Else
            Dim result = New ANVector(Of U, T)(m.Columns)
            For j = 1 To m.Columns
                For i = 1 To m.Rows
                    result(j) = result(j) + v(i) * m(i, j) ''' OPERATOR OVERLOADS BEING USED
                Next
            Next
            Return result
        End If
    End Operator
    'Matrix times Matrix
    Public Shared Operator *(ByVal left As AMNMatrix(Of U, T), ByVal right As AMNMatrix(Of U, T)) _
    As AMNMatrix(Of U, T)
        If left.Columns <> right.Rows Then
            Throw New VectorSpaceException("The number of columns in the left matrix must equal" & _
                                           "The number of rows in the right matrix")
        Else
            Dim result = New AMNMatrix(Of U, T)(left.Rows, right.Columns)
            Dim inner = left.Columns
            For i = 1 To result.Rows
                For j = 1 To result.Columns
                    For k = 1 To inner ''' OPERATOR OVERLOADS BEING USED
                        result(i, j) = result(i, j) + left(i, k) * right(k, j)
                    Next
                Next
            Next
            Return result
        End If
    End Operator
...
End Class

There we go: operator overloading at three levels! I’ve uploaded a zipped Visual-Studio 2008 project to my Public folder with all this code.
Play around, let me know what you think.