Linear Algebra
Colin Fahey
1. Software
LinearAlgebra.zip
Linear algebra source code (C#)
19910 bytes
MD5: 11d8c8035cac30ba543e5e0b72ee9767
2. Introduction
This article describes vectors and matrices in (d)-dimensional space.
3. (d)-dimensional space : attributes
3.1 Array
"array" : A collection of variables such that each variable has a unique name, and such that the names can be assigned an order.
Integer values can be used as names for the variables in an array.
For example, if an array contains (d) variables, then the integer values { 0, 1, 2, ..., (d-1) } can be the names assigned to the variables in the array.
3.2 (d)-dimensional vector
"(d)-dimensional vector" : An array of (d) variables.
"vector component" : A variable within a vector.
3.3 (d)-dimensional vector space
"one-dimensional space" : The complete set of values that can be stored by one variable.
"(d)-dimensional space" : The complete set of combinations of values that can be stored by an array of (d) variables.
The formal definition of a "vector space":
Let (T) be a basic type (examples: real number, integer, complex number, rational number, etc).
Any variable of a basic type is named a "scalar".
A "(T)-type (d)-dimensional vector space" is the set (S) of (d)-dimensional vectors having two operations, vector addition (+), and scalar multiplication (*), satisfying the conditions below:
(1) If (v) and (w) are any two vectors in (S), then (v + w) is also a vector in (S);
(2) If (u), (v), and (w) are any three vectors in (S), then (u + v) + w = u + (v + w);
[additive commutativity]
(3) If (v) and (w) are any two vectors in (S), then (v + w) = (w + v);
[additive associativity]
(4) There is a "zero vector", (0), in (S), such that for any vector (v) in (S), (v + (0)) = v;
[additive identity]
(5) If (c) is any scalar of type (T), and (v) is any vector in (S), then the product (c * v) is a vector in (S);
(6) If (a), (b), and (c) are any scalars of type (T), and (v) and (w) are any vectors in (S), then (a + b) * v = a*v + b*v, and c*(v + w) = c*v + c*w;
[multiplicative distributivity]
(7) If (a) and (b) are any scalars of type (T), and (v) is any vector in (S), then (a*b)*v = a*(b*v);
(8) If "1" is a scalar of type (T) such that (1*1)=1, and (v) is any vector in (S), then (1*v) = v;
(9) For each vector (v) in (S), the vector (-1)*v = -v satisfies v + (-v) = (0);
[additive inverse]
3.4 (d)-dimensional vector code
The code below shows how a (d)-dimensional vector, with 64-bit floating-point components, can be implemented.
An array is sufficient to represent a vector.
The code below has an array contained in a class, only for convenience.
The code is not intended to be efficient.
A structure (example: "struct", a value type) representing vectors of a fixed number of dimensions (examples: 3 or 4) is likely to be much more efficient than the general (d)-dimensional vector class shown here.
Although the code below defines a vector with floating-point components, this document also makes use of vectors with integer components.
The code below can easily be modified to implement vectors with integer components.
using System;
using System.Collections.Generic;
using System.Text;
// Multidimensional vector with 64-bit floating-point components:
public class VectorF64
{
private double[] components;
public int Dimensions()
{
if (null == this.components)
{
return (0);
}
return (this.components.Length);
}
public VectorF64()
{
this.components = null;
}
public VectorF64( params double[] paramValues )
{
this.components = null;
if (null == paramValues)
{
return;
}
int dimensions = paramValues.Length;
this.components = new double[dimensions];
for (int i = 0; i < dimensions; i++)
{
this.components[i] = paramValues[i];
}
}
public VectorF64( VectorF64 other )
{
this.components = null;
if (null == other)
{
return;
}
if (null == other.components)
{
return;
}
int dimensions = other.Dimensions();
this.components = new double[dimensions];
for (int i = 0; i < dimensions; i++)
{
this.components[i] = other.components[i];
}
}
public double this[int index]
{
get
{
if (null == this.components)
{
return (0.0);
}
if
(
(index >= 0)
&& (index < this.components.Length)
)
{
return (this.components[index]);
}
return (0.0);
}
set
{
if (null == this.components)
{
return;
}
if
(
(index >= 0)
&& (index < this.components.Length)
)
{
this.components[index] = value;
}
}
}
public void Write( int precision )
{
if (null == this.components)
{
Log.Write( String.Empty + '(' + ' ' + ')' );
return;
}
int dimensions = this.Dimensions();
if (0 == dimensions)
{
Log.Write( String.Empty + '(' + ' ' + ')' );
return;
}
// Determine the largest component width in characters
// so that we can make all components an equal width.
int largestComponentWidth = 1;
for (int i = 0; i < dimensions; i++)
{
// { index [,minwidth] [:typeCode[precision]] }
// (minwidth<0) means left-justify
String text = String.Format( String.Empty + '{' + '0' + ':'
+ 'g' + precision + '}', this[i] );
if (text.Length > largestComponentWidth)
{
largestComponentWidth = text.Length;
}
}
Log.Write( '(' );
for (int i = 0; i < dimensions; i++)
{
Log.Write( ' ' );
String text =
String.Format( String.Empty + '{' + '0' + ','
+ largestComponentWidth + ':' + 'g' + precision + '}',
this[i] );
Log.Write( text );
if ((i + 1) < dimensions)
{
Log.Write( ',' );
}
else
{
Log.Write( ' ' );
}
}
Log.Write( ')' );
}
public void WriteLine( int precision )
{
this.Write( precision );
Log.WriteLine();
}
public void WriteLine()
{
const int defaultPrecision = 8;
this.Write( defaultPrecision );
Log.WriteLine();
}
// . . .
public static void Test()
{
// A 3-dimensional vector with 64-bit floating-point components:
VectorF64 v3 = new VectorF64( 0.0, 1.0, 2.0 );
v3.WriteLine(); // ( 0, 1, 2 )
// A 4-dimensional vector with 64-bit floating-point components:
VectorF64 v4 = new VectorF64( 0.0, 1.0, 2.0, 3.0 );
v4.WriteLine(); // ( 0, 1, 2, 3 )
// . . .
}
}
"(d)-dimensional zero vector" : A vector with all (d) components equal to zero.
public class VectorF64
{
// . . .
public static VectorF64 Zero( int dimensions )
{
VectorF64 zero = new VectorF64();
zero.components = new double[dimensions];
for (int i = 0; i < dimensions; i++)
{
zero[i] = 0.0;
}
return (zero);
}
// . . .
public static void Test()
{
// . . .
// An 8-dimensional vector with all 8 64-bit floating-point
// components set to zero:
VectorF64 z = VectorF64.Zero( 8 );
z.WriteLine(); // ( 0, 0, 0, 0, 0, 0, 0, 0 )
// . . .
}
}
3.5 (d)-dimensional space : points
"(d)-dimensional space point" : An array of (d) variables with specific values ("coordinate values").
"(d)-dimensional space origin" : An array of (d) variables with all values equal to zero.
3.6 (d)-dimensional vectors : non-relative and relative
A (d)-dimensional vector that is "non-relative" is a (d)-dimensional vector that directly represents a status or configuration.
A (d)-dimensional vector that is "relative" is a (d)-dimensional vector that represents changes to a set of components.
A relative vector can represent the difference between two non-relative vectors.
Given a relative vector, determining a non-relative status or configuration using that relative vector requires combining that relative vector with a non-relative vector.
Non-relative vectors and relative vectors are both vectors.
Whether a particular vector is non-relative or relative must be specified when the vector is defined.
If a (d)-dimensional vector is interpreted as being non-relative, then the (d)-dimensional vector can represent a point in (d)-dimensional space.
3.7 (d)-dimensional vector addition, subtraction, and scaling
Vector addition, subtraction, and scaling:
public class VectorF64
{
// . . .
public static VectorF64 operator +( VectorF64 a, VectorF64 b )
{
if ((null == a) || (null == b))
{
return (new VectorF64()); // Vector not specified.
}
if ((null == a.components) || (null == b.components))
{
return (new VectorF64()); // Vector is empty.
}
if (a.Dimensions() != b.Dimensions())
{
return (new VectorF64()); // Vectors not the same size.
}
int dimensions = a.Dimensions();
VectorF64 result = VectorF64.Zero( dimensions );
for (int i = 0; i < dimensions; i++)
{
result[i] = a[i] + b[i];
}
return (result);
}
public static VectorF64 operator -( VectorF64 a, VectorF64 b )
{
if ((null == a) || (null == b))
{
return (new VectorF64()); // Vector not specified.
}
if ((null == a.components) || (null == b.components))
{
return (new VectorF64()); // Vector is empty.
}
if (a.Dimensions() != b.Dimensions())
{
return (new VectorF64()); // Vectors not the same size.
}
int dimensions = a.Dimensions();
VectorF64 result = VectorF64.Zero( dimensions );
for (int i = 0; i < dimensions; i++)
{
result[i] = a[i] - b[i];
}
return (result);
}
public static VectorF64 operator -( VectorF64 a )
{
if (null == a)
{
return (new VectorF64()); // Vector not specified.
}
if (null == a.components)
{
return (new VectorF64()); // Vector is empty.
}
int dimensions = a.Dimensions();
VectorF64 result = VectorF64.Zero( dimensions );
for (int i = 0; i < dimensions; i++)
{
result[i] = (-( a[i] ));
}
return (result);
}
public static VectorF64 operator *( double scale, VectorF64 a )
{
if (null == a)
{
return (new VectorF64()); // Vector not specified.
}
if (null == a.components)
{
return (new VectorF64()); // Vector is empty.
}
int dimensions = a.Dimensions();
VectorF64 result = VectorF64.Zero( dimensions );
for (int i = 0; i < dimensions; i++)
{
result[i] = scale * a[i];
}
return (result);
}
public static VectorF64 operator *( VectorF64 a, double scale )
{
if (null == a)
{
return (new VectorF64()); // Vector not specified.
}
if (null == a.components)
{
return (new VectorF64()); // Vector is empty.
}
int dimensions = a.Dimensions();
VectorF64 result = VectorF64.Zero( dimensions );
for (int i = 0; i < dimensions; i++)
{
result[i] = scale * a[i];
}
return (result);
}
// . . .
public static void Test()
{
// . . .
// Examples of vector addition, subtraction, and scaling:
VectorF64 a = new VectorF64( 0.0, 1.0, 2.0, 3.0 );
a.WriteLine(); // ( 0, 1, 2, 3 )
VectorF64 b = new VectorF64( 3.0, 2.0, 1.0, 0.0 );
b.WriteLine(); // ( 3, 2, 1, 0 )
VectorF64 c = new VectorF64();
c.WriteLine(); // ( )
c = a + b;
c.WriteLine(); // ( 3, 3, 3, 3 )
c = a - b;
c.WriteLine(); // ( -3, -1, 1, 3 )
c = -b;
c.WriteLine(); // ( -3, -2, -1, 0 )
c = 3.0 * a;
c.WriteLine(); // ( 0, 3, 6, 9 )
// . . .
}
}
3.8 (d)-dimensional basis vectors
The formal definition of a "basis" of a vector space:
Let (T) be a basic type (examples: real number, integer, complex number, rational number, etc).
Any variable of a basic type is named a "scalar".
Let (V) be a "(T)-type (d)-dimensional vector space".
If the non-zero vectors { u1, u2, ..., ud } in (V) are such that every vector (v) in (V) can be written as a "linear combination" of those vectors, v = c1*u1 + c2*u2 + ... + cd*ud, where { c1, c2, ..., cd } are scalars of type (T), then (V) is "spanned" by vectors { u1, u2, ..., ud }.
Any set of non-zero vectors { u1, u2, ..., ud } that "span" vector space (V) is named a "basis" of (V).
One simple "basis" of a (d)-dimensional vector space is the set of (d) distinct (d)-dimensional vectors, each with one component equal to one and all other components equal to zero.
Such basis vectors are "orthonormal", meaning that they are mutually-perpendicular ("orthogonal") and that each vector has unit length.
Each such vector is a unit vector parallel to one of the (d) coordinate axes.
Expressing an arbitrary vector as a "linear combination" of these basis vectors is direct; each component of the arbitrary vector is multiplied by the corresponding basis vector, and these products are added together to form the arbitrary vector.
The code below shows how a vector can be expressed as a "linear combination" of basis vectors.
The code below arbitrarily defines an orthonormal set of basis vectors.
public class VectorF64
{
// . . .
public static VectorF64 BasisVector( int dimensions, int componentIndex )
{
if (dimensions < 0)
{
// Invalid number of dimensions specified.
return (new VectorF64());
}
VectorF64 basisVector = VectorF64.Zero( dimensions );
if ((componentIndex >= 0) && (componentIndex < dimensions))
{
basisVector[ componentIndex ] = 1.0;
}
return (basisVector);
}
// . . .
public static void Test()
{
// . . .
// The 4 basis vectors of 4-dimensional space,
// each with 4 64-bit floating-point components:
VectorF64 b0 = VectorF64.BasisVector( 4, 0 );
b0.WriteLine(); // ( 1, 0, 0, 0 )
VectorF64 b1 = VectorF64.BasisVector( 4, 1 );
b1.WriteLine(); // ( 0, 1, 0, 0 )
VectorF64 b2 = VectorF64.BasisVector( 4, 2 );
b2.WriteLine(); // ( 0, 0, 1, 0 )
VectorF64 b3 = VectorF64.BasisVector( 4, 3 );
b3.WriteLine(); // ( 0, 0, 0, 1 )
// . . .
}
}
Any vector in (d)-dimensional space can be expressed as a sum of the products of numbers and basis vectors:
public class VectorF64
{
// . . .
public static void Test()
{
// . . .
// The following two vectors are equivalent:
// A 4-dimensional vector with 64-bit floating-point components:
VectorF64 va = new VectorF64( 0.1, 1.1, 2.2, 3.3 );
va.WriteLine(); // ( 0.1, 1.1, 2.2, 3.3 )
// A 4-dimensional vector formed by scaling the 4 independent
// basis vectors for 4-dimensional space:
VectorF64 vb =
0.1 * VectorF64.BasisVector( 4, 0 )
+ 1.1 * VectorF64.BasisVector( 4, 1 )
+ 2.2 * VectorF64.BasisVector( 4, 2 )
+ 3.3 * VectorF64.BasisVector( 4, 3 );
vb.WriteLine(); // ( 0.1, 1.1, 2.2, 3.3 )
// . . .
}
}
3.9 (d)-dimensional space : distance between points
Let (P) be a (d)-dimensional vector that represents a point in (d)-dimensional space.
Let (Q) be a (d)-dimensional vector that represents a point in (d)-dimensional space.
Let (R) be a (d)-dimensional vector that represents the change in (d) coordinates to get from point (P) to point (Q); R = (Q - P).
In one-dimensional space, P = ( p0 ), Q = ( q0 ), and R = (Q - P) = ( q0 - p0 ).
The distance between the two points is: Abs( q0 - p0 ).
In two-dimensional space, P = ( p0, p1 ), Q = ( q0, q1 ), and R = (Q - P) = ( q0-p0, q1-p1 ).
Interpreting the two perpendicular displacements as the perpendicular sides of a right-triangle, the distance between the points corresponds to the length of the hypotenuse of that triangle.
The Pythagorean formula, (a*a) + (b*b) = (c*c), where (a) and (b) are the lengths of the perpendicular sides of a right-triangle, and (c) is the length of the hypotenuse (skew side), can be used to determine the distance between the two points: Sqrt( Sq(q0-p0) + Sq(q1-p1) ).
In three-dimensional space, P = ( p0, p1, p2 ), Q = ( q0, q1, q2 ), and R = (Q - P) = ( q0-p0, q1-p1, q2-p2 ).
Interpreting the displacement ( q0-p0, q1-p1, 0 ) as the perpendicular sides of a right-triangle, and using the Pythagorean formula, gives the distance between point (P) and the point ( q0, q1, p2 ): d01 = Sqrt( Sq(q0-p0) + Sq(q1-p1) ).
The displacement ( q0-p0, q1-p1, 0 ) is perpendicular to the displacement ( 0, 0, q2-p2 ), and another right-triangle can be formed, and the Pythagorean formula can be used again.
Thus, the distance from point (P) to point (Q) is given by: Sqrt( Sq(d01) + Sq(q2-p2) ) = Sqrt( (Sq(q0-p0) + Sq(q1-p1)) + Sq(q2-p2) ) = Sqrt( Sq(q0-p0) + Sq(q1-p1) + Sq(q2-p2) ).
The method of extending the distance formula from two-dimensional space to three-dimensional space can be applied repeatedly to eventually determine the distance formula for (d)-dimensional space: Sqrt( (Sq(q0-p0) + Sq(q1-p1)) + Sq(q2-p2) + ... + Sq(qd-pd) ).
The code below defines a function named "Length" which computes the length of a (d)-dimensional vector.
When a vector represents the displacement between two points in (d)-dimensional space, the length of the vector represents the distance between those two points.
public class VectorF64
{
// . . .
public double Length()
{
if (null == this.components)
{
return (0.0); // Vector empty.
}
int dimensions = this.Dimensions();
double sumOfSquares = 0.0;
for (int i = 0; i < dimensions; i++)
{
sumOfSquares += (this[i] * this[i]);
}
double length = Math.Sqrt( sumOfSquares );
return (length);
}
public static double Length( VectorF64 a )
{
if (null == a)
{
return (0.0); // Vector not specified.
}
if (null == a.components)
{
return (0.0); // Vector is empty.
}
return (a.Length());
}
// . . .
public static void Test()
{
// . . .
// Example of vector length:
// A 6-dimensional vector representing a point (p):
VectorF64 p = new VectorF64( 0.0, 1.0, 2.0, 3.0, 4.0, 5.0 );
p.WriteLine(); // ( 0, 1, 2, 3, 4, 5 )
// A 6-dimensional vector representing a point (q):
VectorF64 q = new VectorF64( -5.0, 4.0, -3.0, 2.0, -1.0, 0.0 );
q.WriteLine(); // ( -5, 4, -3, 2, -1, 0 )
// A 6-dimensional vector representing the displacement
// from point (p) to point (q):
VectorF64 r = q - p;
r.WriteLine(); // ( -5, 3, -5, -1, -5, -5 )
// The distance between point (p) and point (q)
// in 6-dimensional space:
double distance = r.Length();
Log.WriteLine( distance ); // 10.4880884817015
// . . .
}
}
3.10 (d)-dimensional vectors : dot product
The "dot product" converts two (d)-dimensional vectors to a number.
The code below computes a dot product of two vectors:
public class VectorF64
{
// . . .
public static double Dot( VectorF64 a, VectorF64 b )
{
if ((null == a) || (null == b))
{
return (0.0); // Vector not specified.
}
if ((null == a.components) || (null == b.components))
{
return (0.0); // Vector is empty.
}
if (a.Dimensions() != b.Dimensions())
{
return (0.0); // Vectors not the same size.
}
int dimensions = a.Dimensions();
double dotProduct = 0.0;
for (int i = 0; i < dimensions; i++)
{
dotProduct += (a[i] * b[i]);
}
return (dotProduct);
}
// . . .
}
For any vector, (A), Length(A) = Sqrt(Dot(A,A)).
3.11 (d)-dimensional vectors : definition of "parallel"
Vectors (A) and (B) are "parallel" if all of the following are true:
(1) A.Length() > 0;
(2) B.Length() > 0;
(3) Abs(Dot(A,B)) = A.Length()*B.Length().
The code below determines if a pair of vectors are parallel (possibly anti-aligned).
Floating-point numbers can accumulate errors in the fractional part due to the limited precision, and, therefore, computer code should include non-zero tolerances when comparing floating-point numbers.
The code includes example tolerance values, but the example tolerance values might not be appropriate for some tasks.
public class VectorF64
{
// . . .
public static bool Parallel( VectorF64 a, VectorF64 b )
{
// Smallest normalized float : 1.1754943e-38
// Smallest normalized double: 2.2250738585072020e-308
double nonZeroThreshold = 1.0e-38; // conservative for double
// double: (52+1)-bit mantissa; log10(2^53)=15.95 decimal digits
double fractionalDifferenceThreshold = 1.0e-14; // conservative
if ((null == a) || (null == b))
{
return (false); // Vector is not specified.
}
if ((null == a.components) || (null == b.components))
{
return (false); // Vector is empty.
}
if (a.Dimensions() != b.Dimensions())
{
return (false); // Vectors not the same size.
}
double lengthA = a.Length();
if (lengthA <= nonZeroThreshold)
{
return (false);
}
double lengthB = b.Length();
if (lengthB <= nonZeroThreshold)
{
return (false);
}
double oneImpliesParallel =
Math.Abs( VectorF64.Dot( a, b ) ) / (lengthA * lengthB);
double absoluteDifferenceFromOne =
Math.Abs( oneImpliesParallel - 1.0 );
if (absoluteDifferenceFromOne <= fractionalDifferenceThreshold)
{
return (true);
}
return (false);
}
// . . .
public static void Test()
{
// . . .
// Example of testing for parallel vectors:
// A 6-dimensional vector:
VectorF64 vf = new VectorF64( 0.0, 1.0, 2.0, 3.0, 4.0, 5.0 );
vf.WriteLine(); // ( 0, 1, 2, 3, 4, 5 )
// A 6-dimensional vector:
VectorF64 vg = new VectorF64( 0.0, -2.0, -4.0, -6.0, -8.0, -10.0 );
vg.WriteLine(); // ( 0, -2, -4, -6, -8, -10 )
// Determine if the specified vectors are parallel
// (or "anti-aligned"):
bool parallel = VectorF64.Parallel( vf, vg );
Log.WriteLine( parallel ); // True
// Add a non-negligible displacement to a component of a vector:
vf[0] += 1.0e-5;
vf.WriteLine(); // ( 1E-05, 1, 2, 3, 4, 5 )
// Determine if the specified vectors are parallel
// (or "anti-aligned"):
parallel = VectorF64.Parallel( vf, vg );
Log.WriteLine( parallel ); // False
// . . .
}
}
3.12 (d)-dimensional vectors : definition of "perpendicular"
Vectors (A) and (B) are "perpendicular" if all of the following are true:
(1) A.Length() > 0;
(2) B.Length() > 0;
(3) Abs(Dot(A,B)) = 0.
The code below determines if a pair of vectors are perpendicular. Floating-point numbers can accumulate errors in the fractional part due to the limited precision, and, therefore, computer code should include non-zero tolerances when comparing floating-point numbers.
The code includes example tolerance values, but the example tolerance values might not be appropriate for some tasks.
public class VectorF64
{
// . . .
public static bool Perpendicular( VectorF64 a, VectorF64 b )
{
// Smallest normalized float : 1.1754943e-38
// Smallest normalized double: 2.2250738585072020e-308
double nonZeroThreshold = 1.0e-38; // conservative for double
// double: (52+1)-bit mantissa; log10(2^53)=15.95 decimal digits
double fractionalDifferenceThreshold = 1.0e-14; // conservative
if ((null == a) || (null == b))
{
return (false); // Vector is not specified.
}
if ((null == a.components) || (null == b.components))
{
return (false); // Vector is empty.
}
if (a.Dimensions() != b.Dimensions())
{
return (false); // Vectors not the same size.
}
double lengthA = a.Length();
if (lengthA <= nonZeroThreshold)
{
return (false);
}
double lengthB = b.Length();
if (lengthB <= nonZeroThreshold)
{
return (false);
}
double zeroImpliesPerpendicular =
Math.Abs( VectorF64.Dot( a, b ) ) / (lengthA * lengthB);
double absoluteDifferenceFromZero =
Math.Abs( zeroImpliesPerpendicular );
if (absoluteDifferenceFromZero <= fractionalDifferenceThreshold)
{
return (true);
}
return (false);
}
// . . .
public static void Test()
{
// . . .
// Example of testing for perpendicular vectors:
// A 6-dimensional vector:
VectorF64 vf2 = new VectorF64( 0.0, 1.0, 2.0, 0.0, 4.0, 5.0 );
vf2.WriteLine(); // ( 0, 1, 2, 0, 4, 5 )
// A 6-dimensional vector:
VectorF64 vg2 = new VectorF64( 10.0, 0.0, 0.0, -5.0, 0.0, 0.0 );
vg2.WriteLine(); // ( 10, 0, 0, -5, 0, 0 )
// Determine if the specified vectors are perpendicular
bool perpendicular = VectorF64.Perpendicular( vf2, vg2 );
Log.WriteLine( perpendicular ); // True
// Add a non-negligible displacement to a component of a vector:
vf2[0] += 1.0e-13;
vf2.WriteLine(); // ( 1E-13, 1, 2, 0, 4, 5 )
// Determine if the specified vectors are perpendicular
perpendicular = VectorF64.Perpendicular( vf2, vg2 );
Log.WriteLine( perpendicular ); // False
// . . .
}
}
3.13 Matrices
"matrix" : A collection of variables such that each variable has a unique combination of a "row" name and a "column" name.
"entry" : A variable within a matrix.
Integer values can be used as "row" names and "column" names for the variables in a matrix.
For example, if a matrix has (totalRows) rows and (totalColumns) columns, then the integer values { 0, 1, ..., (totalRows-1) } can be the names assigned to the rows, and the integer values { 0, 1, ..., (totalColumns-1) } can be the names assigned to the columns.
Thus, a variables in a matrix can be specified by specifying a pair of integers, ( row, column ), indicating the combination of row and column indices corresponding to the specific variable.
The size of a matrix is specified as "(totalRows) * (totalColumns)" (or "(totalRows) by (totalColumns)").
This order of dimensions is the same as the order of dimensions used to specify entries in the matrix ("( row, column )").
[ This convention is somewhat unfortunate, because for many two-dimensional uses (examples: images, graphs, etc) the usual convention is to specify dimensions as "width * height" and coordinates as "( horizontal, vertical )" (or "( x, y )").
This is the opposite of the order of dimensions and coordinates used to describe matrices and their entries. ]
A matrix with (totalRows) equal to (totalColumns) is named "square"; otherwise, the matrix is named "rectangular".
A matrix can be regarded as containing a set of "row vectors", where the variables in each row are interpreted as belonging to a vector.
A matrix can also be regarded as containing a set of "column vectors", where the variables in each column are interpreted as belonging to a vector.
Matrices can represent a wide variety of mathematical relationships.
The meaning of a matrix, and the operations that might be appropriate for processing entries of a matrix, depends on the context.
However, there are basic rules of matrix arithmetic that are relevant to many contexts, and these basic rules will be defined in a subsequent section.
An array and a (totalColumns) value are sufficient to represent a matrix.
The array can have (totalRows * totalColumns) variables, and the entry at ( row, column ) can correspond to the array variable at index ((totalColumns * row) + column).
The code below defines a matrix, with 64-bit floating-point entries.
An array and a (totalColumns) value are sufficient to represent a matrix.
The code below has an array and a (totalColumns) value contained in a class, only for convenience.
The code below is not intended to be efficient.
A structure (example: "struct", a value type) representing matrices of particular fixed dimensions (examples: 2*2, 3*3, or 4*4) is likely to be much more efficient than the general (totalRows * totalColumns) class shown here.
Although the code below defines a matrix with floating-point entries, this article also makes use of matrices with integer entries.
The code below can be easily modified to implement matrices with integer entries.
using System;
using System.Collections.Generic;
using System.Text;
// Matrix with 64-bit floating-point entries:
public class MatrixF64
{
private int totalRows;
private int totalColumns;
private double[] entries;
// matrix[ row, column ] = entries[ (totalColumns * row) + column ]
public int Columns()
{
return (this.totalColumns);
}
public int Rows()
{
return (this.totalRows);
}
public MatrixF64()
{
this.totalRows = 0;
this.totalColumns = 0;
this.entries = null;
}
public MatrixF64( int rows, int columns, params double[] paramValues )
{
this.totalRows = 0;
this.totalColumns = 0;
this.entries = null;
if (rows <= 0)
{
return; // Specified an invalid row count.
}
if (columns <= 0)
{
return; // Specified an invalid column count.
}
int totalEntries = (rows * columns);
this.totalRows = rows;
this.totalColumns = columns;
this.entries = new double[totalEntries];
for (int k = 0; k < totalEntries; k++)
{
this.entries[k] = 0.0;
}
if (null == paramValues)
{
return; // No entries specified.
}
int totalParamValues = paramValues.Length;
if (totalParamValues != totalEntries)
{
// The number of specified entries does not match the
// total number of matrix entries.
return;
}
for (int k = 0; k < totalParamValues; k++)
{
this.entries[k] = paramValues[k];
}
}
public MatrixF64( MatrixF64 other )
{
this.totalRows = 0;
this.totalColumns = 0;
this.entries = null;
if (null == other)
{
return;
}
if
(
(null == other.entries)
|| (0 == other.totalRows)
|| (0 == other.totalColumns)
)
{
return;
}
int totalEntries = (other.totalRows * other.totalColumns);
if (other.entries.Length != totalEntries)
{
return;
}
this.totalRows = other.totalRows;
this.totalColumns = other.totalColumns;
this.entries = new double[totalEntries];
for (int k = 0; k < totalEntries; k++)
{
this.entries[k] = other.entries[k];
}
}
public double this[int row, int column]
{
get
{
if
(
(this.totalRows > 0)
&& (this.totalColumns > 0)
&& (null != this.entries)
&& (row >= 0)
&& (row < this.totalRows)
&& (column >= 0)
&& (column < this.totalColumns)
)
{
int k = (this.totalColumns * row) + column;
if ((k >= 0) && (k < this.entries.Length))
{
return (this.entries[k]);
}
}
return (0.0);
}
set
{
if
(
(this.totalRows > 0)
&& (this.totalColumns > 0)
&& (null != this.entries)
&& (row >= 0)
&& (row < this.totalRows)
&& (column >= 0)
&& (column < this.totalColumns)
)
{
int k = (this.totalColumns * row) + column;
if ((k >= 0) && (k < this.entries.Length))
{
this.entries[k] = value;
}
}
}
}
public void WriteLine( int precision )
{
if
(
(this.totalRows <= 0)
|| (this.totalColumns <= 0)
|| (null == this.entries)
)
{
Log.WriteLine( String.Empty + '[' + ' ' + ']' );
return;
}
// Determine the largest entry width in characters
// so that we can make all columns an equal width.
int largestEntryWidth = 1;
for (int i = 0; i < this.totalRows; i++)
{
for (int j = 0; j < this.totalColumns; j++)
{
// { index [,minwidth] [:typeCode[precision]] }
// (minwidth<0) means left-justify
String text = String.Format( String.Empty
+ '{' + '0' + ':' + 'g' + precision + '}', this[i,j] );
if (text.Length > largestEntryWidth)
{
largestEntryWidth = text.Length;
}
}
}
// Print each row of the matrix.
for (int i = 0; i < this.totalRows; i++)
{
Log.Write( '[' );
for (int j = 0; j < this.totalColumns; j++)
{
Log.Write( ' ' );
String text = String.Format( String.Empty
+ '{' + '0' + ',' + largestEntryWidth + ':' + 'g'
+ precision + '}', this[i,j] );
Log.Write( text );
if ((j + 1) < this.totalColumns)
{
Log.Write( ',' );
}
else
{
Log.Write( ' ' );
}
}
Log.WriteLine( ']' );
}
}
public void WriteLine()
{
const int defaultPrecision = 8;
this.WriteLine( defaultPrecision );
}
// . . .
public static void Test()
{
// A 5x3 matrix (5 rows x 3 columns) with
// 64-bit floating-point entries:
MatrixF64 m =
new MatrixF64
(
5, 3,
0.0, 1.0, 2.0, // row 0
3.0, 4.0, 5.0, // row 1
6.0, 7.0, 8.0, // row 2
9.0, 10.0, 11.0, // row 3
12.0, 13.0, 14.0 // row 4
);
m.WriteLine();
// [ 0, 1, 2 ] // row 0
// [ 3, 4, 5 ] // row 1
// [ 6, 7, 8 ] // row 2
// [ 9, 10, 11 ] // row 3
// [ 12, 13, 14 ] // row 4
// . . .
}
}
"zero matrix" : A matrix with all entries equal to zero.
public class MatrixF64
{
// . . .
public static MatrixF64 Zero( int rows, int columns )
{
MatrixF64 zero = new MatrixF64( rows, columns );
// The constructor above sets all entries to zero.
return (zero);
}
// . . .
public static void Test()
{
// . . .
// An 8x2 matrix (8 rows x 2 columns) with all 16
// 64-bit floating-point entries set to zero:
MatrixF64 z = MatrixF64.Zero( 8, 2 );
z.WriteLine();
// [ 0, 0 ] // row 0
// [ 0, 0 ] // row 1
// [ 0, 0 ] // row 2
// [ 0, 0 ] // row 3
// [ 0, 0 ] // row 4
// [ 0, 0 ] // row 5
// [ 0, 0 ] // row 6
// [ 0, 0 ] // row 7
// . . .
}
}
3.14 Matrix addition, subtraction, and multiplication
Matrix addition, subtraction, and multiplication.
public class MatrixF64
{
// . . .
public static MatrixF64 operator +( MatrixF64 a, MatrixF64 b )
{
if ((null == a) || (null == b))
{
return (null);
}
MatrixF64 result = MatrixF64.Zero( a.totalRows, a.totalColumns );
if
(
(null == a.entries)
|| (null == b.entries)
|| (a.totalRows != b.totalRows)
|| (a.totalColumns != b.totalColumns)
)
{
return (result);
}
for (int i = 0; i < a.totalRows; i++)
{
for (int j = 0; j < a.totalColumns; j++)
{
result[i, j] = a[i, j] + b[i, j];
}
}
return (result);
}
public static MatrixF64 operator -( MatrixF64 a, MatrixF64 b )
{
if ((null == a) || (null == b))
{
return (null);
}
MatrixF64 result = MatrixF64.Zero( a.totalRows, a.totalColumns );
if
(
(null == a.entries)
|| (null == b.entries)
|| (a.totalRows != b.totalRows)
|| (a.totalColumns != b.totalColumns)
)
{
return (result);
}
for (int i = 0; i < a.totalRows; i++)
{
for (int j = 0; j < a.totalColumns; j++)
{
result[i, j] = a[i, j] - b[i, j];
}
}
return (result);
}
public static MatrixF64 operator -( MatrixF64 a )
{
if (null == a)
{
return (null);
}
MatrixF64 result = MatrixF64.Zero( a.totalRows, a.totalColumns );
if (null == a.entries)
{
return (result);
}
for (int i = 0; i < a.totalRows; i++)
{
for (int j = 0; j < a.totalColumns; j++)
{
result[i, j] = (-( a[i, j] ));
}
}
return (result);
}
public static MatrixF64 operator *( double scale, MatrixF64 a )
{
if (null == a)
{
return (null);
}
MatrixF64 result = MatrixF64.Zero( a.totalRows, a.totalColumns );
if (null == a.entries)
{
return (result);
}
for (int i = 0; i < a.totalRows; i++)
{
for (int j = 0; j < a.totalColumns; j++)
{
result[i, j] = scale * a[i, j];
}
}
return (result);
}
public static MatrixF64 operator *( MatrixF64 a, double scale )
{
if (null == a)
{
return (null);
}
MatrixF64 result = MatrixF64.Zero( a.totalRows, a.totalColumns );
if (null == a.entries)
{
return (result);
}
for (int i = 0; i < a.totalRows; i++)
{
for (int j = 0; j < a.totalColumns; j++)
{
result[i, j] = scale * a[i, j];
}
}
return (result);
}
public static MatrixF64 operator *( MatrixF64 a, MatrixF64 b )
{
if ((null == a) || (null == b))
{
return (null);
}
MatrixF64 result = MatrixF64.Zero( a.totalRows, b.totalColumns );
// The number of columns of the first matrix must be equal
// to the number of rows of the second matrix.
if
(
(null == a.entries)
|| (null == b.entries)
|| (a.totalColumns != b.totalRows)
)
{
return (result);
}
// Entry (i,j) of the result is equal to the dot product
// of row (i) of (a) and column (j) of (b).
for (int i = 0; i < a.totalRows; i++)
{
for (int j = 0; j < b.totalColumns; j++)
{
double dotProduct = 0.0;
for (int k = 0; k < a.totalColumns; k++)
{
dotProduct += (a[i, k] * b[k, j]);
}
result[i, j] = dotProduct;
}
}
return (result);
}
// . . .
public static void Test()
{
// . . .
// Examples of matrix addition, subtraction, and multiplication:
MatrixF64 a =
new MatrixF64
(
3, 5,
0.0, 1.0, 2.0, 3.0, 4.0, // row 0
5.0, 6.0, 7.0, 8.0, 9.0, // row 1
10.0, 11.0, 12.0, 13.0, 14.0 // row 2
);
MatrixF64 b =
new MatrixF64
(
3, 5,
4.0, -3.0, 2.0, -1.0, 0.0, // row 0
-9.0, 8.0, -7.0, 6.0, -5.0, // row 1
14.0, -13.0, 12.0, -11.0, 10.0 // row 2
);
MatrixF64 c =
new MatrixF64
(
5, 3,
0.0, 1.0, 2.0, // row 0
3.0, 4.0, 5.0, // row 1
6.0, 7.0, 8.0, // row 2
9.0, 10.0, 11.0, // row 3
12.0, 13.0, 14.0 // row 4
);
MatrixF64 result = new MatrixF64();
result = a + b;
result.WriteLine();
// [ 4, -2, 4, 2, 4 ]
// [ -4, 14, 0, 14, 4 ]
// [ 24, -2, 24, 2, 24 ]
result = a - b;
result.WriteLine();
// [ -4, 4, 0, 4, 4 ]
// [ 14, -2, 14, 2, 14 ]
// [ -4, 24, 0, 24, 4 ]
result = -a;
result.WriteLine();
// [ 0, -1, -2, -3, -4 ]
// [ -5, -6, -7, -8, -9 ]
// [ -10, -11, -12, -13, -14 ]
result = 3.0 * a;
result.WriteLine();
// [ 0, 3, 6, 9, 12 ]
// [ 15, 18, 21, 24, 27 ]
// [ 30, 33, 36, 39, 42 ]
result = a * c; // (3x5) * (5x3) = (3x3)
result.WriteLine();
// [ 90, 100, 110 ]
// [ 240, 275, 310 ]
// [ 390, 450, 510 ]
result = c * a; // (5x3) * (3x5) = (5x5)
result.WriteLine();
// [ 25, 28, 31, 34, 37 ]
// [ 70, 82, 94, 106, 118 ]
// [ 115, 136, 157, 178, 199 ]
// [ 160, 190, 220, 250, 280 ]
// [ 205, 244, 283, 322, 361 ]
// . . .
}
}
"identity matrix" : A square matrix (total rows equals total columns) with entries equal to (1) on the diagonal (entries with a row index equal to its column index), and with all other entries equal to (0).
If a square matrix, (M), is multiplied by an "identity matrix", (I), of the same number of rows (and columns), the product is equal to (M).
Multiplying (I) by (M) also produces a product equal to (M).
Thus, the "identity matrix" is analogous to the number "1" for the multiplication of numbers (scalars).
The code below creates an identity matrix with a specified number of rows.
public class MatrixF64
{
// . . .
public static MatrixF64 Identity( int rows )
{
MatrixF64 identity = MatrixF64.Zero( rows, rows );
for (int i = 0; i < rows; i++)
{
identity[i, i] = 1.0;
}
return (identity);
}
// . . .
public static void Test()
{
// . . .
// Examples of multiplying by an identity matrix:
MatrixF64 identity3x3 = MatrixF64.Identity( 3 );
// [ 1, 0, 0 ]
// [ 0, 1, 0 ]
// [ 0, 0, 1 ]
MatrixF64 s3x3 =
new MatrixF64
(
3, 3,
0.0, 1.0, 2.0, // row 0
3.0, 4.0, 5.0, // row 1
6.0, 7.0, 8.0 // row 2
);
result = s3x3 * identity3x3;
result.WriteLine();
// [ 0, 1, 2 ]
// [ 3, 4, 5 ]
// [ 6, 7, 8 ]
// (the same as s3x3)
result = identity3x3 * s3x3;
result.WriteLine();
// [ 0, 1, 2 ]
// [ 3, 4, 5 ]
// [ 6, 7, 8 ]
// (the same as s3x3)
// . . .
}
}
3.15 Matrix : LU factoring; back-substitution
"Matrix LU factoring" : A procedure that converts a square matrix, (M), to two new square matrices, (L) and (U), having the same size as (M), such that (L) * (U) = (M), and such that matrix (U) is "upper triangular" (all entries below the diagonal are zero), and matrix (L) is "lower triangular" (all entries above the diagonal are zero).
Matrix LU factoring can be used as part of a larger procedure to solve a system of equations, or to find the inverse of a matrix, or to find the determinant of a matrix.
Matrix LU factoring is an alternative to the "Gauss-Jordan elimination" procedure.
Gauss-Jordan elimination requires a system of equations (A)*(x)=(b), whereas LU factoring only requires a matrix (A).
Also, after determining the LU factoring of a matrix (A), it is very easy to determine (x) given any (b).
The procedure to solve for the vector (x) in (A)*(x)=(b), given (b) and the LU factoring (L)*(U)=(A), involves a procedure that will be named "back-substitution" here.
The code below includes a procedure for computing the LU factors of any square, non-singular matrix.
The code below includes a procedure for doing LU back-substitution.
Caution: The computer code below computes the L and U factors of a row-permuted version of a specified matrix.
Also, the L and U matrix results are combined in to a single output matrix.
public class MatrixF64
{
// . . .
public static void FindLUFactorsOfARowPermutedVersionOfAnOriginalMatrix
(
MatrixF64 originalMatrix,
ref MatrixF64 LUFactorsMatrix,
ref int[] rowPermutationOfOriginalMatrix,
ref bool rowExchangeParityIsOdd
)
{
// ''LU factoring'' involves factoring a square matrix (M) in to a
// lower-triangular matrix (L), and an upper-triangular matrix (U),
// such that (L)*(U) = (M).
//
// However, the specified matrix (originalMatrix) might not be
// optimal for direct LU factoring, due to the locations of extreme
// values in the matrix. Therefore, this function instead implicitly
// factors a row-permuted version of originalMatrix, where the
// permutation of rows is determined by the values in the specified
// matrix. The specific row permutation selected by this function is
// returned in an array of row indices (rowPermutationOfOriginalMatrix).
// The resulting (L) and (U) factors are merged in to a single
// output matrix (LUFactorsMatrix). Therefore:
//
// (L part of LUFactorsMatrix) * (U part of LUFactorsMatrix)
// = (originalMatrix with rows permuted according to
// rowPermutationOfOriginalMatrix).
//
// Although factoring a row-permuted version of originalMatrix
// makes it more difficult to interpret the results of this function,
// the resulting LU factoring is likely to be more accurate. In a
// sense, this function indicates (via rowPermutationOfOriginalMatrix)
// how to permute the rows of originalMatrix to be able to produce
// the most accurate LU factoring directly, and this function produces
// that factoring (via the LUFactorsMatrix output).
//
// When using the (L) and (U) matrices (merged in to LUFactorsMatrix)
// to solve a system of equations, it is necessary to permute the rows
// of the solution vector of the system of equations according to
// rowPermutation.
//
// The output matrix (LUFactorsMatrix) is a combination of
// the (L) and (U) matrices:
//
// To form the (L) matrix from (LUFactorsMatrix):
// assume (0) for entries above the diagonal,
// assume (1) for entries on the diagonal,
// and use (LUFactorsMatrix) entries below the diagonal.
//
// To form the (U) matrix from (LUFactorsMatrix):
// assume (0) for entries below the diagonal,
// and use (LUFactorsMatrix) entries on and above the diagonal.
//
// The output array (rowPermutationOfOriginalMatrix) indicates how the
// rows of the original matrix have implicitly (not actually) been
// permuted.
//
// The output boolean (rowExchangeParityIsOdd) indicates if the parity
// of the permutation of rows is odd.
//
// This implementation is partly based on pseudo-code appearing in
// ''Introduction to Algorithms'' (Cormen; Lieserson; Rivest;
// 24th printing, 2000; MIT Press; ISBN 0-262-03141-8).
// See the section named ''Overview of LUP decomposition'', in
// chapter 31 (''Matrix Operations''). The implementation here follows
// the ''LUP-Decomposition'' pseudo-code appearing in a section named
// ''Computing an LUP decomposition''.
//
// This implementation also uses ideas found in the implementation of
// ''ludcmp'' in the book ''Numerical recipes in C : The art of
// scientific computing'' (Press; Teukolsky; Vetterling; Flannery;
// second edition; 1996 reprinting; Cambridge University Press).
// The overall loop structure here resembles that of ''ludcmp''.
// The idea of determining and caching row scale factors in advance of
// finding pivots is used here. The method in ''ludcmp'' to avoid zero
// pivot values has been modified here to handle excessively-small
// pivot values, too. HOWEVER, the row permutation indices produced by
// the following implementation is a true permutation of
// {0,1,...,(totalRows-1)}, whereas ''ludcmp'' (and ''lubksb'')
// interpret their ''permutation indices'' as swaps to perform
// (roughly, for each (i), swap b[ i ] and b[ indx[i] ]). So, the
// following implementation of LU-factoring is not directly compatible
// with ''ludcmp'' or ''lubksb''. Also, the following procedure
// does not destroy the original matrix (unlike ''ludcmp'').
double singularMatrixIfMaxRowElementIsLessThanThis = (1.0e-19);
double forcePivotsToHaveAtLeastThisAbsoluteValue = (1.0e-19);
LUFactorsMatrix = null;
rowPermutationOfOriginalMatrix = null;
rowExchangeParityIsOdd = false;
if (null == originalMatrix)
{
return; // No matrix specified.
}
if
(
(originalMatrix.totalRows <= 0)
|| (originalMatrix.totalColumns <= 0)
|| (null == originalMatrix.entries)
)
{
return; // Matrix is empty.
}
if (originalMatrix.totalRows != originalMatrix.totalColumns)
{
return; // Matrix is not square.
}
// Duplicate the original matrix
LUFactorsMatrix = new MatrixF64( originalMatrix );
// (Lines 2-3 of LUP-Decomposition)
// Initialize the row permutation array to the identity
// permutation.
rowPermutationOfOriginalMatrix = new int[LUFactorsMatrix.totalRows];
for (int i = 0; i < LUFactorsMatrix.totalRows; i++)
{
rowPermutationOfOriginalMatrix[i] = i;
}
// For each row, determine the largest absolute value of
// the row elements, and use the reciprocal as the scale for
// that row.
VectorF64 rowScales = VectorF64.Zero( LUFactorsMatrix.totalRows );
for (int i = 0; i < LUFactorsMatrix.totalRows; i++)
{
double largestElementAbsoluteValueInRow = 0.0;
for (int j = 0; j < LUFactorsMatrix.totalColumns; j++)
{
double absoluteElementValue =
Math.Abs( LUFactorsMatrix[i, j] );
if (absoluteElementValue >
largestElementAbsoluteValueInRow)
{
largestElementAbsoluteValueInRow =
absoluteElementValue;
}
}
if (largestElementAbsoluteValueInRow <
singularMatrixIfMaxRowElementIsLessThanThis)
{
return; // Matrix is singular
}
rowScales[i] = (1.0 / largestElementAbsoluteValueInRow);
}
// (Lines 4-18 of LUP-Decomposition)
// Go through the columns of the matrix.
for (int j = 0; j < LUFactorsMatrix.totalColumns; j++)
{
// (Lines 17-18 of LUP-Decomposition)
// (or equation (2.3.12) of NR)
for (int i = 0; i < j; i++)
{
double sum = LUFactorsMatrix[i, j];
for (int k = 0; k < i; k++)
{
sum -= (LUFactorsMatrix[i, k] * LUFactorsMatrix[k, j]);
}
LUFactorsMatrix[i, j] = sum;
}
// Go through all remaining rows and find the best pivot.
double largestScaledSum = 0.0;
int rowIndexOfLargestScaledSum = j;
for (int i = j; i < LUFactorsMatrix.totalRows; i++)
{
// (equation (2.3.13) of NR)
double sum = LUFactorsMatrix[i, j];
for (int k = 0; k < j; k++)
{
sum -= (LUFactorsMatrix[i, k] * LUFactorsMatrix[k, j]);
}
LUFactorsMatrix[i, j] = sum;
double scaledSum = rowScales[i] * Math.Abs( sum );
if (scaledSum >= largestScaledSum)
{
largestScaledSum = scaledSum;
rowIndexOfLargestScaledSum = i;
}
}
// If indeed we found a better pivot, then exchange rows.
if (j != rowIndexOfLargestScaledSum)
{
// (Line 12 of LUP-Decomposition)
// Exchange the row permutation indices
int tempRowIndex = rowPermutationOfOriginalMatrix[j];
rowPermutationOfOriginalMatrix[j] =
rowPermutationOfOriginalMatrix[rowIndexOfLargestScaledSum];
rowPermutationOfOriginalMatrix[rowIndexOfLargestScaledSum] =
tempRowIndex;
// (Lines 13-14 of LUP-Decomposition)
// Exchange the elements of the rows
for (int k = 0; k < LUFactorsMatrix.totalColumns; k++)
{
double temp = LUFactorsMatrix[rowIndexOfLargestScaledSum, k];
LUFactorsMatrix[rowIndexOfLargestScaledSum, k] =
LUFactorsMatrix[j, k];
LUFactorsMatrix[j, k] = temp;
}
// Exchange the row scale factors
double scaleFactor =
rowScales[rowIndexOfLargestScaledSum];
rowScales[rowIndexOfLargestScaledSum] = rowScales[j];
rowScales[j] = scaleFactor;
// Invert the overall row exchange parity
rowExchangeParityIsOdd = (!(rowExchangeParityIsOdd));
}
// Force the pivot element to have at least a certain
// absolute value.
if (Math.Abs( LUFactorsMatrix[j, j] ) <
forcePivotsToHaveAtLeastThisAbsoluteValue)
{
if (LUFactorsMatrix[j, j] < 0.0)
{
LUFactorsMatrix[j, j] =
(-(forcePivotsToHaveAtLeastThisAbsoluteValue));
}
else
{
LUFactorsMatrix[j, j] =
forcePivotsToHaveAtLeastThisAbsoluteValue;
}
}
// If not the final column, then divide all column elements
// below the diagonal by the pivot element, matrixLU[j,j].
// (Lines 15-16 of LUP-Decomposition)
if (j != (LUFactorsMatrix.totalColumns - 1))
{
double reciprocalOfPivot = (1.0 / LUFactorsMatrix[j, j]);
for (int i = (j + 1); i < LUFactorsMatrix.totalRows; i++)
{
LUFactorsMatrix[i, j] *= reciprocalOfPivot;
}
}
}
}
public static void LUBacksubstitution
(
MatrixF64 LUFactorsMatrix,
int[] rowPermutationOfOriginalMatrix,
VectorF64 givenProductVector,
ref VectorF64 solutionVector
)
{
// This implementation is based on pseudo-code appearing in
// ''Introduction to Algorithms'' (Cormen; Lieserson; Rivest;
// 24th printing, 2000; MIT Press; ISBN 0-262-03141-8).
// See the section named ''Overview of LUP decomposition'', in
// chapter 31 (''Matrix Operations''). The implementation here
// follows the ''LUP-Solve'' pseudo-code appearing in a
// section named ''forward and back substitution''.
solutionVector = null;
if (null == LUFactorsMatrix)
{
return; // No matrix specified.
}
if (null == rowPermutationOfOriginalMatrix)
{
return; // No permutation specified.
}
if (null == givenProductVector)
{
return; // No product vector specified.
}
if
(
(null == LUFactorsMatrix.entries)
|| (LUFactorsMatrix.totalRows <= 0)
|| (LUFactorsMatrix.totalColumns <= 0)
)
{
return; // Matrix is empty.
}
if (LUFactorsMatrix.totalRows !=
LUFactorsMatrix.totalColumns)
{
return; // Matrix is not square.
}
if (givenProductVector.Dimensions() !=
LUFactorsMatrix.totalRows)
{
return; // Product vector size not equal to matrix row count.
}
// Copy the product vector in to the result.
solutionVector = new VectorF64( givenProductVector );
// (See LUP-Solve, lines 2-3)
for (int i = 0; i < LUFactorsMatrix.totalRows; i++)
{
double sum =
givenProductVector[rowPermutationOfOriginalMatrix[i]];
for (int j = 0; j < i; j++)
{
sum -= (LUFactorsMatrix[i, j] * solutionVector[j]);
}
solutionVector[i] = sum;
}
// (See LUP-Solve, lines 4-5)
for (int i = (LUFactorsMatrix.totalRows - 1); i >= 0; i--)
{
double sum = solutionVector[i];
for
(
int j = (i + 1);
j < LUFactorsMatrix.totalColumns;
j++
)
{
sum -= (LUFactorsMatrix[i, j] * solutionVector[j]);
}
double diagonalElement = LUFactorsMatrix[i, i];
if (diagonalElement != 0.0)
{
solutionVector[i] = (sum / diagonalElement);
}
}
}
// . . .
public static void Test()
{
// . . .
// Example of LU factoring and back-substitution:
MatrixF64 a3x3 =
new MatrixF64
(
3, 3,
0.0, 1.0, 2.0, // row 0
3.0, 4.0, 5.0, // row 1
6.0, 7.0, 8.0 // row 2
);
int[] rowPermutationOfOriginalMatrix = null;
bool rowInterchangeParityIsOdd = false;
MatrixF64 a3x3LU = null;
MatrixF64.FindLUFactorsOfARowPermutedVersionOfAnOriginalMatrix
(
a3x3,
ref a3x3LU,
ref rowPermutationOfOriginalMatrix,
ref rowInterchangeParityIsOdd
);
a3x3LU.WriteLine();
// [ 6, 7, 8 ]
// [ 0, 1, 2 ]
// [ 0.5, 0.5, 1e-19 ]
// The matrix above is a combination of the L and U
// matrices. The L matrix is 0 above the diagonal,
// 1 on the diagonal, and the shown values below the
// diagonal. The U matrix is 0 below the diagonal,
// and the shown values on and above the diagonal.
// This LU factoring is of a ROW-PERMUTED version
// of the original matrix.
// The following L and U matrices are manually formed
// by inspecting the LU factoring matrix above.
MatrixF64 L3x3 =
new MatrixF64
(
3, 3,
1.0, 0.0, 0.0, // row 0
0.0, 1.0, 0.0, // row 1
0.5, 0.5, 1.0 // row 2
);
MatrixF64 U3x3 =
new MatrixF64
(
3, 3,
6.0, 7.0, 8.0, // row 0
0.0, 1.0, 2.0, // row 1
0.0, 0.0, 0.0 // row 2
);
// The product (L)*(U) produces a ROW-PERMUTED version
// of the original matrix:
result = L3x3 * U3x3;
result.WriteLine( 4 );
// [ 6, 7, 8 ]
// [ 0, 1, 2 ]
// [ 3, 4, 5 ]
VectorF64 b3x1 = new VectorF64( 1.0, 2.0, 3.0 );
VectorF64 x3x1 = null;
MatrixF64.LUBacksubstitution
(
a3x3LU,
rowPermutationOfOriginalMatrix,
b3x1,
ref x3x1
);
// Solution vector
x3x1.WriteLine();
// ( -0.66666667, 1, 0 )
// Both the given b3x1 and this x3x1 are relative to the
// original matrix, a3x3. The a3x3LU for a ROW-PERMUTED
// version of a3x3, but LUBacksubstitution accepts a
// non-permuted product vector (e.g., b3x1) and produces
// a non-permuted solution vector (e.g., x3x1).
// Check solution by multiplying the original matrix a3x3
// by the solution vector x3x1, to get a product vector.
VectorF64 checkb3x1 = a3x3 * x3x1;
checkb3x1.WriteLine();
// ( 1, 2, 3 )
// This matches the b3x1 vector we specified above, so
// the solution is valid.
// . . .
}
}
3.16 Matrix : determinant
"determinant" : A function that converts any square (n * n) matrix (A) to a number, det(A), such that:
(1) If a matrix (B) results from exchanging two rows, or two columns, of a matrix (A), then det(B) = (-(det(A)));
(2) If a matrix (B) results from multiplying any row, or any column, of a matrix (A), by a number (c), then det(B) = c * det(A);
(3) If a matrix (B) results from adding a multiple of one row to another row, or from adding a multiple of one column to another column, then det(B) = det(A).
(4) If (A) is an upper-triangular matrix (A[i,j]=0 for all (i>j)) or a lower-triangular matrix (A[i,j]=0 for all (i<j)), then det(A) = (A[0,0] * A[1,1] * ... * A[n-1,n-1]);
( For example, the determinant of an identity matrix is one; det(I) = 1. )
Rules (1), (2), and (3), can be used, in a process named "Gaussian elimination", to convert any square matrix to a triangular matrix.
Rule (4) can be used to compute the determinant of a triangular matrix.
The formal definition of a "determinant":
Let M[n](K) denote the set of all (n * n) matrices over the field (K).
The "determinant" is the unique function (F) such that F:M[n](K) --> K with two attributes:
(1) F is alternating multilinear with regard to columns (or rows); and,
(2) F( I ) = 1.
A "multilinear" function (D) of an (n * n) matrix (A) can be written as: D(A) = Sum( A[0,k[0]] * A[1,k[1]] * ... * A[n-1,k[n-1]] * D( e[k[0]], e[k[1]], ..., e[k[n-1]] ) ), where the sum is taken over all (n^n) combinations of 0 <= k[i] <= (n-1), and where e[j] represents row (j) of the identity matrix.
Such a function can be characterized by the values of D( e[k[0]], e[k[1]], ..., e[k[n-1]] ).
An "alternating multilinear" function is a multilinear function, (D), such that D( ..., e[i], ..., e[i], ... ) = 0, and D( ..., e[i], e[j], ... ) = (-( D( ..., e[j], e[i], ... ) )).
Thus, swapping columns (or rows) changes the sign of the function.
Also, the term D( e[k[0]], ..., e[k[n-1]] ) is zero if the combination of k[i] values does not involve all unique values; if the set of values is not a permutation of the numbers {0,1,...,(n-1)}.
Letting D( I ) = 1 (D( e[0], e[1], ..., e[n-1] ) = 1), in combination with the alternating multilinear attributes described above, means that all values of the alternating multilinear function D( e[k[0]], e[k[1]], ..., e[k[n-1]] ) can be determined.
The function will be non-zero only if the set of k[i] values is a permutation of {0,1,...,(n-1)}, and will have a magnitude of one if it is non-zero, and the sign will be the number of swaps required to convert the permutation to the identity permutation.
The value of a determinant of an (n * n) matrix (A) can be computed by adding up products of the form (sgn(p) * A[0,p[0]] * A[1,p[1]] * ... * A[(n-1),p[n-1]]) for each and every permutation (p) of the numbers {0,1,2,...,(n-1)}.
The term (sgn(p)) denotes the "signature" (or swap count) of permutation (p), where (sgn(p)) is equal to (+1) if (p) is an "even permutation", and is equal to (-1) if (p) is an "odd permutation".
Because there are (n!) permutations of the numbers {0,1,2,...,(n-1)}, this formula has (n!) summands.
Other procedures for computing the determinant (examples: those involving Gaussian elimination, or LU factoring, or expansion by minors, etc) implicitly compute a hierarchy of intermediate quantities that enable those procedures to compute the determinant in an order of (n^3) steps.
If the determinant of a matrix is zero, the matrix does not have an inverse.
If the determinant of a matrix is non-zero, the matrix has an inverse.
If a matrix represents the coefficients of a linear system of equations, and the determinant of that matrix is zero, then the system of equations does not have a unique solution.
If a matrix represents the coefficients of a linear system of equations, and the determinant of that matrix is non-zero, then the system of equations has a unique solution.
The code below includes procedures for computing the determinant of any square matrix.
Determinants for 1*1, 2*2, 3*3, and 4*4 matrices are computed by explicit formulas.
Using explicit formulas to compute the determinants for small matrices is likely to compute results faster than using a general procedure to compute those determinants.
The procedure for the general case uses LU factoring, but there are many other methods of computing determinants.
public class MatrixF64
{
// . . .
private double Determinant1x1()
{
if
(
(this.totalRows != 1)
|| (this.totalColumns != 1)
|| (null == this.entries)
)
{
return (0.0);
}
// The determinant of a 1x1 matrix is simply
// the value of the single element.
return ( this[0,0] );
}
private double Determinant2x2()
{
if
(
(this.totalRows != 2)
|| (this.totalColumns != 2)
|| (null == this.entries)
)
{
return (0.0);
}
return
(
+ this[0,0] * this[1,1]
- this[0,1] * this[1,0]
);
}
private double Determinant3x3()
{
if
(
(this.totalRows != 3)
|| (this.totalColumns != 3)
|| (null == this.entries)
)
{
return (0.0);
}
return
(
+ this[0,0] * this[1,1] * this[2,2]
+ this[0,1] * this[1,2] * this[2,0]
+ this[0,2] * this[1,0] * this[2,1]
- this[0,0] * this[1,2] * this[2,1]
- this[0,1] * this[1,0] * this[2,2]
- this[0,2] * this[1,1] * this[2,0]
);
}
private double Determinant4x4()
{
if
(
(this.totalRows != 4)
|| (this.totalColumns != 4)
|| (null == this.entries)
)
{
return (0.0);
}
return
(
+ this[0, 0] * this[1, 1] * this[2, 2] * this[3, 3]
+ this[0, 0] * this[1, 2] * this[2, 3] * this[3, 1]
+ this[0, 0] * this[1, 3] * this[2, 1] * this[3, 2]
+ this[0, 1] * this[1, 0] * this[2, 3] * this[3, 2]
+ this[0, 1] * this[1, 2] * this[2, 0] * this[3, 3]
+ this[0, 1] * this[1, 3] * this[2, 2] * this[3, 0]
+ this[0, 2] * this[1, 0] * this[2, 1] * this[3, 3]
+ this[0, 2] * this[1, 1] * this[2, 3] * this[3, 0]
+ this[0, 2] * this[1, 3] * this[2, 0] * this[3, 1]
+ this[0, 3] * this[1, 0] * this[2, 2] * this[3, 1]
+ this[0, 3] * this[1, 1] * this[2, 0] * this[3, 2]
+ this[0, 3] * this[1, 2] * this[2, 1] * this[3, 0]
- this[0, 0] * this[1, 1] * this[2, 3] * this[3, 2]
- this[0, 0] * this[1, 2] * this[2, 1] * this[3, 3]
- this[0, 0] * this[1, 3] * this[2, 2] * this[3, 1]
- this[0, 1] * this[1, 0] * this[2, 2] * this[3, 3]
- this[0, 1] * this[1, 2] * this[2, 3] * this[3, 0]
- this[0, 1] * this[1, 3] * this[2, 0] * this[3, 2]
- this[0, 2] * this[1, 0] * this[2, 3] * this[3, 1]
- this[0, 2] * this[1, 1] * this[2, 0] * this[3, 3]
- this[0, 2] * this[1, 3] * this[2, 1] * this[3, 0]
- this[0, 3] * this[1, 0] * this[2, 1] * this[3, 2]
- this[0, 3] * this[1, 1] * this[2, 2] * this[3, 0]
- this[0, 3] * this[1, 2] * this[2, 0] * this[3, 1]
);
}
public double Determinant()
{
if
(
(this.totalRows <= 0)
|| (this.totalColumns <= 0)
|| (null == this.entries)
)
{
return (0.0); // Matrix is empty.
}
if (this.totalRows != this.totalColumns)
{
// Matrix is not square
return (0.0);
}
// Simple cases
if (1 == this.totalRows) { return (this.Determinant1x1()); }
if (2 == this.totalRows) { return (this.Determinant2x2()); }
if (3 == this.totalRows) { return (this.Determinant3x3()); }
if (4 == this.totalRows) { return (this.Determinant4x4()); }
int[] rowPermutationOfOriginalMatrix = null;
bool rowExchangeParityIsOdd = false;
MatrixF64 originalMatrix = new MatrixF64( this );
MatrixF64 LUFactorsMatrix = null;
MatrixF64.FindLUFactorsOfARowPermutedVersionOfAnOriginalMatrix
(
originalMatrix,
ref LUFactorsMatrix,
ref rowPermutationOfOriginalMatrix,
ref rowExchangeParityIsOdd
);
if ((null == LUFactorsMatrix)
|| (null == rowPermutationOfOriginalMatrix))
{
return (0.0);
}
double determinant = 1.0;
if (true == rowExchangeParityIsOdd)
{
determinant = (-1.0);
}
for (int j = 0; j < LUFactorsMatrix.totalColumns; j++)
{
determinant *= LUFactorsMatrix[j, j];
}
return (determinant);
}
// . . .
public static void Test()
{
// . . .
// Examples of matrix determinants:
MatrixF64 d2x2 =
new MatrixF64
(
2, 2,
2.0, 5.0, // row 0
-4.0, -7.0 // row 1
);
double detd2x2 = d2x2.Determinant();
Log.WriteLine( detd2x2 );
// detd2x2 = 6
MatrixF64 d3x3 =
new MatrixF64
(
3, 3,
2.0, 5.0, 7.0, // row 0
-4.0, -1.0, 6.0, // row 1
9.0, 8.0, 3.0 // row 2
);
double detd3x3 = d3x3.Determinant();
Log.WriteLine( detd3x3 );
// detd3x3 = 67
MatrixF64 d4x4 =
new MatrixF64
(
4, 4,
7.0, -5.0, 2.0, 4.0, // row 0
3.0, 2.0, 6.0, 3.0, // row 1
-9.0, 8.0, -3.0, 2.0, // row 2
5.0, 3.0, 2.0, 5.0 // row 3
);
double detd4x4 = d4x4.Determinant();
Log.WriteLine( detd4x4 );
// detd4x4 = 1457
// . . .
}
}
3.17 (d)-dimensional vectors : cross product
"cross product" : A function associated with a (d)-dimensional space that converts an ordered list of (d-1) (d)-dimensional vectors, { A(0), A(1), ..., A(d-2) }, to a (d)-dimensional vector, such that the function has the attributes below:
(1) If A(i)=0 (zero vector), for any 0 <= i <= (d-2), then Cross( A(0), A(1), ..., A(d-2) )=0 (zero vector).
The cross product of an ordered list of vectors is zero if any vector in that list is zero;
(2) If any pair of vectors, A(i) and A(j), for any (i!=j) with 0 <= i <= (d-2) and 0 <= j <= (d-2), are parallel (Abs(Dot(A(i),A(j))) = Length(A(i)) * Length(A(j))), then Cross( A(0), A(1), ..., A(d-2) )=0 (zero vector).
The cross product of an ordered list of vectors is zero if any two vectors in that list are parallel;
(3) If every vector A(i) is non-zero, and if A(i) is not parallel with A(j) for all (i!=j) with 0 <= i <= (d-2) and 0 <= j <= (d-2), then B = Cross( A(0), A(1), ..., A(d-2) ) is such that Dot(A(i),B)=0 for all 0 <= i <= (d-2).
If the cross product of an ordered list of vectors is non-zero, then the cross product result is perpendicular to every vector in the ordered list of vectors;
(4) If (c) is a numerical value, then Cross( (c * A(0)), A(1), ..., A(d-2) ) = Cross( A(0), (c * A(1)), ..., A(d-2) ) = ... = Cross( A(0), A(1), ..., (c * A(d-2)) ) = c * Cross( A(0), A(1), ..., A(d-2) ).
Given an ordered list of vectors, the cross product of the ordered list of vectors with any one vector multiplied by a specific numerical value is equivalent to multiplying the cross product of the ordered list of vectors by the numerical value;
(5) Given two numbers (i) and (j), where (i!=j) and 0 <= i <= (d-2) and 0 <= j <= (d-2), and an ordered list of (d-1) (d)-dimensional vectors, { A(0), A(1), ..., A(d-2) }, and another ordered list of (d-1) (d)-dimensional vectors, { B(0), B(1), ..., B(d-2) }, such that B(k) = A(k) for all ((k != i) && (k != j)) with 0 <= k <= (d-2), and such that B(i) = A(j) and B(j) = A(i), then Cross( A(0), A(1), ..., A(d-2) ) = (-1) * Cross( B(0), B(1), ..., B(d-2) ).
If the cross product of an ordered list of vectors is non-zero, then that cross product has the opposite sign of the cross product of that same ordered list of vectors with the exception of exactly two vectors being exchanged (swapped).
Given basis vectors e(k), for 0 <= k <= (d-1), of (d)-dimensional space, the cross product of an ordered list of (d-1) vectors can be computed by computing the determinant of the (d*d) matrix below:
e(0), e(1), ..., e(d-1),
A( 0 )[0], A( 0 )[1], ..., A( 0 )[d-1],
A( 1 )[0], A( 1 )[1], ..., A( 1 )[d-1],
...,
A(d-2)[0], A(d-2)[1], ..., A(d-2)[d-1]
The determinant of that matrix will not be a number, but will instead be a vector, where there will be numerical coefficients to the (d) basis vectors.
For two-dimensional space, and the vector A = ( ax, ay ):
Cross( A )
= Determinant
(
e(x), e(y),
ax, ay
)
= ay * e(x)
- ax * e(y)
= ( ay, -ax )
The result vector is perpendicular to the single input vector.
If the x-y plane is viewed such that x increases in a rightward direction and y increases in an upward direction, then the cross product result is a quarter-turn "clockwise" relative to the input vector.
Examples of cross products in two-dimensional space:
Cross( e(x) ) = (-1) * e(y)
Cross( e(y) ) = e(x)
For three-dimensional space, and the vectors A = ( ax, ay, az ) and B = ( bx, by, bz ):
Cross( A, B )
= Determinant
(
e(x), e(y), e(z),
ax, ay, az,
bx, by, bz
)
= (ay*bz - az*by) * e(x)
+ (az*bx - ax*bz) * e(y)
+ (ax*by - ay*bx) * e(z)
= ( (ay*bz - az*by), (az*bx - ax*bz), (ax*by - ay*bx) )
Examples of cross products in three-dimensional space:
Cross( e(x), e(y) ) = e(z)
Cross( e(y), e(x) ) = (-1) * e(z)
Cross( e(y), e(z) ) = e(x)
Cross( e(z), e(y) ) = (-1) * e(x)
Cross( e(z), e(x) ) = e(y)
Cross( e(x), e(z) ) = (-1) * e(y)
For four-dimensional space, and the vectors A = ( ax, ay, az, aw ), B = ( bx, by, bz, bw ), and C = ( cx, cy, cz, cw ):
Cross( A, B, C )
= Determinant
(
e(x), e(y), e(z), e(w),
ax, ay, az, aw,
bx, by, bz, bw,
cx, cy, cz, cw
)
= ( ay*bz*cw - ay*bw*cz + az*bw*cy - az*by*cw + aw*by*cz - aw*bz*cy) * e(x)
+ (-ax*bz*cw + ax*bw*cz - az*bw*cx + az*bx*cw - aw*bx*cz + aw*bz*cx) * e(y)
+ ( ax*by*cw - ax*bw*cy + ay*bw*cx - ay*bx*cw + aw*bx*cy - aw*by*cx) * e(z)
+ (-ax*by*cz + ax*bz*cy - ay*bz*cx + ay*bx*cz - az*bx*cy + az*by*cx) * e(w)
= ( ( ay*bz*cw - ay*bw*cz + az*bw*cy - az*by*cw + aw*by*cz - aw*bz*cy),
(-ax*bz*cw + ax*bw*cz - az*bw*cx + az*bx*cw - aw*bx*cz + aw*bz*cx),
( ax*by*cw - ax*bw*cy + ay*bw*cx - ay*bx*cw + aw*bx*cy - aw*by*cx),
(-ax*by*cz + ax*bz*cy - ay*bz*cx + ay*bx*cz - az*bx*cy + az*by*cx) )
Examples of cross products in four-dimensional space:
Cross( e(x), e(y), e(z) ) = (-1) * e(w)
Cross( e(x), e(z), e(y) ) = e(w)
Cross( e(y), e(x), e(z) ) = e(w)
Cross( e(y), e(z), e(x) ) = (-1) * e(w)
Cross( e(z), e(x), e(y) ) = (-1) * e(w)
Cross( e(z), e(y), e(x) ) = e(w)
Cross( e(x), e(y), e(w) ) = e(z)
Cross( e(x), e(w), e(y) ) = (-1) * e(z)
Cross( e(y), e(x), e(w) ) = (-1) * e(z)
Cross( e(y), e(w), e(x) ) = e(z)
Cross( e(w), e(x), e(y) ) = e(z)
Cross( e(w), e(y), e(x) ) = (-1) * e(z)
Cross( e(x), e(z), e(w) ) = (-1) * e(y)
Cross( e(x), e(w), e(z) ) = e(y)
Cross( e(z), e(x), e(w) ) = e(y)
Cross( e(z), e(w), e(x) ) = (-1) * e(y)
Cross( e(w), e(x), e(z) ) = (-1) * e(y)
Cross( e(w), e(z), e(x) ) = e(y)
Cross( e(y), e(z), e(w) ) = e(x)
Cross( e(y), e(w), e(z) ) = (-1) * e(x)
Cross( e(z), e(y), e(w) ) = (-1) * e(x)
Cross( e(z), e(w), e(y) ) = e(x)
Cross( e(w), e(y), e(z) ) = e(x)
Cross( e(w), e(z), e(y) ) = (-1) * e(x)
The code below computes the cross-product of two or more vectors.
The test examples show that the non-zero result vectors are perpendicular to all of the input vectors.
public class MatrixF64
{
// . . .
public static VectorF64 Cross( params VectorF64[] parameterArray )
{
if (null == parameterArray)
{
return(new VectorF64());
}
int totalVectors = parameterArray.Length;
// The number of dimensions of the space must be one larger
// than the number of specified vectors.
int dimensions = (1 + totalVectors);
if (dimensions < 2)
{
// The cross product does not exist for space with
// fewer than two dimensions.
return (new VectorF64());
}
// All specified vectors must exist and must have a number
// of components equal to the number of dimensions of the space.
for (int k = 0; k < totalVectors; k++)
{
if (null == parameterArray[k])
{
// Vector missing.
return (new VectorF64());
}
else if (null == parameterArray[k].components)
{
// Vector is empty.
return (new VectorF64());
}
else if (parameterArray[k].components.Length != dimensions)
{
// Vector has the wrong number of components.
return (new VectorF64());
}
}
// Form a (d*d) matrix with the (d-1) supplied vectors
// forming the final (d-1) rows of the matrix.
MatrixF64 matrix = MatrixF64.Zero( dimensions, dimensions );
for (int i = 1; i < dimensions; i++)
{
VectorF64 rowVector = parameterArray[ i - 1 ];
for (int j = 0; j < dimensions; j++)
{
matrix[i, j] = rowVector[j];
}
}
VectorF64 result = VectorF64.Zero( dimensions );
// For each component of the result vector, set the first
// row of the matrix to a basis vector, compute the matrix
// determinant, and use the resulting value as the
// component of the result vector. This is inefficient,
// due to the determinant being computed (d) times.
// Clearly it would be better to have explicit formulas
// for 2-dimensional, 3-dimensional, and 4-dimensional
// cross-products.
for (int i = 0; i < dimensions; i++)
{
// Set the first row of the matrix to be the basis
// vector for the coordinate axis with index (i).
for (int j = 0; j < dimensions; j++)
{
matrix[0, j] = 0.0;
}
matrix[0, i] = 1.0;
// Compute the determinant of the modified matrix.
double determinant = matrix.Determinant();
// Set component (i) of the result to the determinant
// that was computed.
result[i] = determinant;
}
return (result);
}
// . . .
public static void Test()
{
// . . .
// Examples of cross products:
// Three-dimensional example:
VectorF64 vcp1x3a = new VectorF64( 1.0, 2.0, 3.0 );
VectorF64 vcp1x3b = new VectorF64( 3.0, 1.0, 2.0 );
VectorF64 vcp1x3 = VectorF64.Cross( vcp1x3a, vcp1x3b );
vcp1x3.WriteLine(); // ( 1, 7, -5 )
// Verify that the result is perpendicular to the original
// vectors:
Log.WriteLine( VectorF64.Perpendicular( vcp1x3, vcp1x3a ) );
Log.WriteLine( VectorF64.Perpendicular( vcp1x3, vcp1x3b ) );
// True, True
// Four-dimensional example:
VectorF64 vcp1x4a = new VectorF64( 1.0, 2.0, 3.0, 4.0 );
VectorF64 vcp1x4b = new VectorF64( 4.0, 1.0, 2.0, 3.0 );
VectorF64 vcp1x4c = new VectorF64( 3.0, 4.0, 1.0, 2.0 );
VectorF64 vcp1x4 = VectorF64.Cross( vcp1x4a, vcp1x4b, vcp1x4c );
vcp1x4.WriteLine(); // ( 4, 4, 44, -36 )
// Verify that the result is perpendicular to the original
// vectors:
Log.WriteLine( VectorF64.Perpendicular( vcp1x4, vcp1x4a ) );
Log.WriteLine( VectorF64.Perpendicular( vcp1x4, vcp1x4b ) );
Log.WriteLine( VectorF64.Perpendicular( vcp1x4, vcp1x4c ) );
// True, True, True
// . . .
}
}
3.18 Matrix : inverse
"matrix inverse" : A square matrix related to another square matrix of the same size such that the product of the two matrices is the "identity matrix".
If a matrix has an associated inverse matrix, there is only one such inverse matrix.
If the determinant of a matrix is zero, then the matrix is named "singular", and the matrix does not have an associated inverse matrix.
The code below computes the inverse of any square, non-singular matrix.
Inverses for 1*1, 2*2, 3*3, and 4*4 matrices are computed by explicit formulas.
Using explicit formulas might be faster than using a general procedure.
The explicit formulas are also reliable.
The procedure for the general case uses LU factoring, but there are many other methods of computing inverses of matrices.
public class MatrixF64
{
// . . .
private MatrixF64 Inverse1x1()
{
if
(
(this.totalRows != 1)
|| (this.totalColumns != 1)
|| (null == this.entries)
)
{
return (MatrixF64.Zero( 1, 1 )); // Matrix is empty.
}
if (0.0 == this[0,0])
{
return ( MatrixF64.Zero( 1, 1 ) ); // Matrix has no inverse.
}
MatrixF64 inverse = MatrixF64.Zero( 1, 1 );
// The inverse of a 1x1 matrix is simply the
// reciprocal of the single entry.
inverse[0,0] = (1.0 / this[0,0]);
return (inverse);
}
private MatrixF64 Inverse2x2()
{
if
(
(this.totalRows != 2)
|| (this.totalColumns != 2)
|| (null == this.entries)
)
{
return (MatrixF64.Zero( 2, 2 )); // Matrix is empty.
}
double determinant = this.Determinant2x2();
if (0.0 == determinant)
{
return (MatrixF64.Zero( 2, 2 )); // Matrix has no inverse.
}
MatrixF64 inverse = MatrixF64.Zero( 2, 2 );
inverse[0,0] = this[1,1];
inverse[0,1] = (-(this[0,1]));
inverse[1,0] = (-(this[1,0]));
inverse[1,1] = this[0,0];
double factor = (1.0 / determinant);
for ( int i = 0; i < inverse.totalRows; i++ )
{
for ( int j = 0; j < inverse.totalColumns; j++ )
{
inverse[i,j] *= factor;
}
}
return (inverse);
}
private MatrixF64 Inverse3x3()
{
if
(
(this.totalRows != 3)
|| (this.totalColumns != 3)
|| (null == this.entries)
)
{
return (MatrixF64.Zero( 3, 3 )); // Matrix is empty.
}
double determinant = this.Determinant3x3();
if (0.0 == determinant)
{
return (MatrixF64.Zero( 3, 3 )); // Matrix has no inverse.
}
MatrixF64 inverse = MatrixF64.Zero( 3, 3 );
inverse[0, 0] =
this[1, 1] * this[2, 2]
- this[1, 2] * this[2, 1];
inverse[0, 1] =
this[0, 2] * this[2, 1]
- this[0, 1] * this[2, 2];
inverse[0, 2] =
this[0, 1] * this[1, 2]
- this[0, 2] * this[1, 1];
inverse[1, 0] =
this[1, 2] * this[2, 0]
- this[1, 0] * this[2, 2];
inverse[1, 1] =
this[0, 0] * this[2, 2]
- this[0, 2] * this[2, 0];
inverse[1, 2] =
this[0, 2] * this[1, 0]
- this[0, 0] * this[1, 2];
inverse[2, 0] =
this[1, 0] * this[2, 1]
- this[1, 1] * this[2, 0];
inverse[2, 1] =
this[0, 1] * this[2, 0]
- this[0, 0] * this[2, 1];
inverse[2, 2] =
this[0, 0] * this[1, 1]
- this[0, 1] * this[1, 0];
double factor = (1.0 / determinant);
for (int i = 0; i < inverse.totalRows; i++)
{
for (int j = 0; j < inverse.totalColumns; j++)
{
inverse[i, j] *= factor;
}
}
return (inverse);
}
private MatrixF64 Inverse4x4()
{
if
(
(this.totalRows != 4)
|| (this.totalColumns != 4)
|| (null == this.entries)
)
{
return (MatrixF64.Zero( 4, 4 )); // Matrix is empty.
}
double determinant = this.Determinant4x4();
if (0.0 == determinant)
{
return (MatrixF64.Zero( 4, 4 )); // Matrix has no inverse.
}
MatrixF64 inverse = MatrixF64.Zero( 4, 4 );
inverse[0,0] =
+ this[1,1] * this[2,2] * this[3,3]
+ this[1,2] * this[2,3] * this[3,1]
+ this[1,3] * this[2,1] * this[3,2]
- this[1,1] * this[2,3] * this[3,2]
- this[1,2] * this[2,1] * this[3,3]
- this[1,3] * this[2,2] * this[3,1];
inverse[0,1] =
+ this[0,1] * this[2,3] * this[3,2]
+ this[0,2] * this[2,1] * this[3,3]
+ this[0,3] * this[2,2] * this[3,1]
- this[0,1] * this[2,2] * this[3,3]
- this[0,2] * this[2,3] * this[3,1]
- this[0,3] * this[2,1] * this[3,2];
inverse[0,2] =
+ this[0,1] * this[1,2] * this[3,3]
+ this[0,2] * this[1,3] * this[3,1]
+ this[0,3] * this[1,1] * this[3,2]
- this[0,1] * this[1,3] * this[3,2]
- this[0,2] * this[1,1] * this[3,3]
- this[0,3] * this[1,2] * this[3,1];
inverse[0,3] =
+ this[0,1] * this[1,3] * this[2,2]
+ this[0,2] * this[1,1] * this[2,3]
+ this[0,3] * this[1,2] * this[2,1]
- this[0,1] * this[1,2] * this[2,3]
- this[0,2] * this[1,3] * this[2,1]
- this[0,3] * this[1,1] * this[2,2];
inverse[1,0] =
+ this[1,0] * this[2,3] * this[3,2]
+ this[1,2] * this[2,0] * this[3,3]
+ this[1,3] * this[2,2] * this[3,0]
- this[1,0] * this[2,2] * this[3,3]
- this[1,2] * this[2,3] * this[3,0]
- this[1,3] * this[2,0] * this[3,2];
inverse[1,1] =
+ this[0,0] * this[2,2] * this[3,3]
+ this[0,2] * this[2,3] * this[3,0]
+ this[0,3] * this[2,0] * this[3,2]
- this[0,0] * this[2,3] * this[3,2]
- this[0,2] * this[2,0] * this[3,3]
- this[0,3] * this[2,2] * this[3,0];
inverse[1,2] =
+ this[0,0] * this[1,3] * this[3,2]
+ this[0,2] * this[1,0] * this[3,3]
+ this[0,3] * this[1,2] * this[3,0]
- this[0,0] * this[1,2] * this[3,3]
- this[0,2] * this[1,3] * this[3,0]
- this[0,3] * this[1,0] * this[3,2];
inverse[1,3] =
+ this[0,0] * this[1,2] * this[2,3]
+ this[0,2] * this[1,3] * this[2,0]
+ this[0,3] * this[1,0] * this[2,2]
- this[0,0] * this[1,3] * this[2,2]
- this[0,2] * this[1,0] * this[2,3]
- this[0,3] * this[1,2] * this[2,0];
inverse[2,0] =
+ this[1,0] * this[2,1] * this[3,3]
+ this[1,1] * this[2,3] * this[3,0]
+ this[1,3] * this[2,0] * this[3,1]
- this[1,0] * this[2,3] * this[3,1]
- this[1,1] * this[2,0] * this[3,3]
- this[1,3] * this[2,1] * this[3,0];
inverse[2,1] =
+ this[0,0] * this[2,3] * this[3,1]
+ this[0,1] * this[2,0] * this[3,3]
+ this[0,3] * this[2,1] * this[3,0]
- this[0,0] * this[2,1] * this[3,3]
- this[0,1] * this[2,3] * this[3,0]
- this[0,3] * this[2,0] * this[3,1];
inverse[2,2] =
+ this[0,0] * this[1,1] * this[3,3]
+ this[0,1] * this[1,3] * this[3,0]
+ this[0,3] * this[1,0] * this[3,1]
- this[0,0] * this[1,3] * this[3,1]
- this[0,1] * this[1,0] * this[3,3]
- this[0,3] * this[1,1] * this[3,0];
inverse[2,3] =
+ this[0,0] * this[1,3] * this[2,1]
+ this[0,1] * this[1,0] * this[2,3]
+ this[0,3] * this[1,1] * this[2,0]
- this[0,0] * this[1,1] * this[2,3]
- this[0,1] * this[1,3] * this[2,0]
- this[0,3] * this[1,0] * this[2,1];
inverse[3,0] =
+ this[1,0] * this[2,2] * this[3,1]
+ this[1,1] * this[2,0] * this[3,2]
+ this[1,2] * this[2,1] * this[3,0]
- this[1,0] * this[2,1] * this[3,2]
- this[1,1] * this[2,2] * this[3,0]
- this[1,2] * this[2,0] * this[3,1];
inverse[3,1] =
+ this[0,0] * this[2,1] * this[3,2]
+ this[0,1] * this[2,2] * this[3,0]
+ this[0,2] * this[2,0] * this[3,1]
- this[0,0] * this[2,2] * this[3,1]
- this[0,1] * this[2,0] * this[3,2]
- this[0,2] * this[2,1] * this[3,0];
inverse[3,2] =
+ this[0,0] * this[1,2] * this[3,1]
+ this[0,1] * this[1,0] * this[3,2]
+ this[0,2] * this[1,1] * this[3,0]
- this[0,0] * this[1,1] * this[3,2]
- this[0,1] * this[1,2] * this[3,0]
- this[0,2] * this[1,0] * this[3,1];
inverse[3,3] =
+ this[0,0] * this[1,1] * this[2,2]
+ this[0,1] * this[1,2] * this[2,0]
+ this[0,2] * this[1,0] * this[2,1]
- this[0,0] * this[1,2] * this[2,1]
- this[0,1] * this[1,0] * this[2,2]
- this[0,2] * this[1,1] * this[2,0];
double factor = (1.0 / determinant);
for (int i = 0; i < inverse.totalRows; i++)
{
for (int j = 0; j < inverse.totalColumns; j++)
{
inverse[i, j] *= factor;
}
}
return (inverse);
}
public MatrixF64 Inverse()
{
if
(
(this.totalRows <= 0)
|| (this.totalColumns <= 0)
|| (null == this.entries)
)
{
return (MatrixF64.Zero( 1, 1 )); // Matrix is empty.
}
if (this.totalRows != this.totalColumns)
{
// Matrix is not square
return (MatrixF64.Zero( 1, 1 ));
}
// Simple cases
if (1 == this.totalRows) { return (this.Inverse1x1()); }
if (2 == this.totalRows) { return (this.Inverse2x2()); }
if (3 == this.totalRows) { return (this.Inverse3x3()); }
if (4 == this.totalRows) { return (this.Inverse4x4()); }
int[] rowPermutationOfOriginalMatrix = null;
bool rowExchangeParityIsOdd = false;
MatrixF64 originalMatrix = new MatrixF64( this );
MatrixF64 LUFactorsMatrix = null;
MatrixF64.FindLUFactorsOfARowPermutedVersionOfAnOriginalMatrix
(
originalMatrix,
ref LUFactorsMatrix,
ref rowPermutationOfOriginalMatrix,
ref rowExchangeParityIsOdd
);
if ((null == LUFactorsMatrix)
|| (null == rowPermutationOfOriginalMatrix))
{
return (MatrixF64.Zero( this.totalRows, this.totalColumns ));
}
VectorF64 columnVector = VectorF64.Zero( totalRows );
VectorF64 solutionVector = VectorF64.Zero( totalRows );
MatrixF64 inverse =
MatrixF64.Zero( this.totalRows, this.totalColumns );
for (int j = 0; j < inverse.totalColumns; j++)
{
for (int i = 0; i < inverse.totalRows; i++)
{
columnVector[i] = 0.0;
}
columnVector[j] = 1.0;
MatrixF64.LUBacksubstitution
(
LUFactorsMatrix,
rowPermutationOfOriginalMatrix,
columnVector,
ref solutionVector
);
if (null != solutionVector)
{
for (int i = 0; i < inverse.totalRows; i++)
{
inverse[i, j] = solutionVector[i];
}
}
}
return (inverse);
}
// . . .
public static void Test()
{
// . . .
// Examples of matrix inverses:
MatrixF64 m2x2 =
new MatrixF64
(
2, 2,
2.0, 5.0, // row 0
-4.0, -7.0 // row 1
);
MatrixF64 inversem2x2 = m2x2.Inverse();
inversem2x2.WriteLine();
// [ -1.1666667, -0.83333333 ]
// [ 0.66666667, 0.33333333 ]
MatrixF64 prodinvm2x2andm2x2 = inversem2x2 * m2x2;
prodinvm2x2andm2x2.WriteLine();
// [ 1, 8.8817842e-16 ]
// [ 0, 1 ]
// This product is very close to a
// 2x2 identity matrix, as it should be.
MatrixF64 prodm2x2andinvm2x2 = m2x2 * inversem2x2;
prodm2x2andinvm2x2.WriteLine();
// [ 1, 0 ]
// [ 0, 1 ]
// This product is the 2x2 identity matrix,
// as it should be.
MatrixF64 m3x3 =
new MatrixF64
(
3, 3,
2.0, 5.0, 7.0, // row 0
-4.0, -1.0, 6.0, // row 1
9.0, 8.0, 3.0 // row 2
);
MatrixF64 inversem3x3 = m3x3.Inverse();
inversem3x3.WriteLine();
// [ -0.76119403, 0.6119403, 0.55223881 ]
// [ 0.98507463, -0.85074627, -0.59701493 ]
// [ -0.34328358, 0.43283582, 0.26865672 ]
MatrixF64 prodinvm3x3andm3x3 = inversem3x3 * m3x3;
prodinvm3x3andm3x3.WriteLine();
// [ 1, 0, 4.4408921e-16 ]
// [ 0, 1, 0 ]
// [ 0, 0, 1 ]
// This product is very close to a 3x3 identity matrix,
// as it should be.
MatrixF64 prodm3x3andinvm3x3 = m3x3 * inversem3x3;
prodm3x3andinvm3x3.WriteLine();
// [ 1, -4.4408921e-16, 4.4408921e-16 ]
// [ 0, 1, -2.220446e-16 ]
// [ -4.4408921e-16, 0, 1 ]
// This product is very close to a 3x3 identity matrix,
// as it should be.
// . . .
}
}
3.19 Multiplicative product of a matrix and a vector
The multiplicative product of a matrix (M) and a vector (a) is a vector (b); (M)*(a) = (b).
Matrix (M) "transforms" vector (a) to vector (b).
The code below computes the multiplicative product of a matrix and a vector.
public class MatrixF64
{
// . . .
public static VectorF64 operator *( MatrixF64 m, VectorF64 v )
{
if (null == m)
{
return (null); // Matrix is not specified.
}
if
(
(null == m.entries)
|| (m.totalRows <= 0)
|| (m.totalColumns <= 0)
)
{
return (null); // Matrix is empty.
}
if (null == v)
{
// Vector is not specified.
return (VectorF64.Zero( m.totalRows ));
}
if (v.Dimensions() != m.totalColumns)
{
// The number of columns of the matrix must be equal
// to the number of rows of the vector.
return (VectorF64.Zero( m.totalRows ));
}
// The result vector has the same number of rows as the matrix.
VectorF64 result = VectorF64.Zero( m.totalRows );
// Component (i) of the result vector is equal to the
// dot product of row (i) of the matrix with the vector.
for (int i = 0; i < m.totalRows; i++)
{
double dotProduct = 0.0;
for (int j = 0; j < m.totalColumns; j++)
{
dotProduct += (m[i, j] * v[j]);
}
result[i] = dotProduct;
}
return (result);
}
// . . .
public static void Test()
{
// . . .
// Examples of a product of a matrix and a vector:
MatrixF64 t3x3 =
new MatrixF64
(
3, 3,
2.0, 5.0, 7.0, // row 0
-4.0, -1.0, 6.0, // row 1
9.0, 8.0, 3.0 // row 2
);
VectorF64 p3x1 = new VectorF64( 1.0, 2.0, 3.0 );
VectorF64 q3x1 = t3x3 * p3x1;
q3x1.WriteLine();
// ( 33, 12, 34 )
// The multiplicative product of the inverse
// of the matrix and the result vector computed
// above should reproduce the original vector.
MatrixF64 inverset3x3 = t3x3.Inverse();
VectorF64 prodinvt3x3andq3x1 = inverset3x3 * q3x1;
prodinvt3x3andq3x1.WriteLine();
// ( 1, 2, 3 )
// This vector is equal to the original vector,
// as expected.
// . . .
}
}
3.20 Using homogeneous matrices to translate homogeneous vectors
"translation" : An operation T(s), specified by a particular vector (s), that converts any vector (r) to a corresponding vector (r + s); T(s) * (r) = (r + s).
T(-s) is the inverse of T(s); Inverse(T(s)) = T(-s); T(-s) * (r) = (r - s).
A translation operation can be represented by a matrix, such that multiplying the translation matrix by a vector has the effect of translating that vector.
The contributions of a translation operation to the components of the result vector do not depend, directly or indirectly, on the components of the original vector.
Thus, a translation matrix must be structured in a manner that enables it to contribute a translation to the result without being affected, directly or indirectly, by the components of the original vector.
A "homogeneous vector" is a vector with a final component equal to 1.
A "homogeneous matrix" is a matrix with a final row with all entries equal to 0 except for a final entry of 1.
The multiplicative product of a homogeneous matrix and another homogeneous matrix is a homogeneous matrix.
The multiplicative product of a homogeneous matrix and a homogeneous vector is a homogeneous vector.
One interesting attribute of this product is that the entries in the final column of the homogeneous matrix, with the exception of the final entry, are added directly to the components of the result vector, without being affected, directly or indirectly, by the components of the original vector.
Thus, a translation operation can be represented by a homogeneous matrix, and such a matrix can be used to translate homogeneous vectors.
If we wish to represent positions and directions in (d)-dimensional space, and we wish to use matrices to translate positions, then homogeneous matrices and homogeneous vectors can be used.
However, the homogeneous vectors must have (d+1) components, and the homogeneous matrices must have ((d+1)*(d+1)) entries.
The final component of a homogeneous vector and the final row of a homogeneous matrix do not correspond to any of the (d) coordinates that we wish to represent in (d)-dimensional space, but instead serve only to make it possible to get matrices to do some operations independently from the (d) coordinate values.
Using homogeneous matrices and homogeneous vectors makes it easier to compound transformations and apply transformations to vectors, but other operations become more inconvenient.
Representing positions or directions in (d)-dimensional space using homogeneous vectors means that the vectors will have (d+1) components, with the final component equal to 1.
Computing the dot product of two homogeneous vectors using a regular dot-product function would combine implicitly the dot product of the first (d) components with an irrelevant value of 1 due to final components of 1.
Therefore, to get the dot product of the relevant (d) components of the homogeneous vector (with (d+1) components), the dot product of homogeneous vectors should be done using a special "homogeneous vector dot product" function that simply skips the final component of the vectors.
Alternatively, the dot product could be computed with the regular dot product function with the understanding that an amount of 1 should be subtracted from the result before interpreting the result as pertaining to the (d) relevant coordinates.
The code below shows how a matrix can be formed that can then be used to translate a vector.
public class MatrixF64
{
// . . .
public static MatrixF64 FormHomogeneousTranslationMatrix( VectorF64 v )
{
// This method assumes the specified vector is a homogeneous
// vector, where the final component is not relevant to the
// translation. Thus, we form a square homogeneous matrix
// with a total number of rows equal to the number of
// components of the specified vector. We form an identity
// matrix of the required size, and copy all but the final
// component of the specified vector to the final column of
// the matrix.
if (null == v)
{
return (MatrixF64.Zero( 1, 1 )); // Vector is not specified.
}
if (v.Dimensions() <= 0)
{
return (MatrixF64.Zero( 1, 1 )); // Vector is empty.
}
MatrixF64 result =
MatrixF64.Identity( v.Dimensions() );
// Add all but the final component of the specified vector
// to the final column of the result matrix. The final
// entry of the final column of the matrix will remain one (1).
for (int i = 0; i < (result.totalRows - 1); i++)
{
result[i, (result.totalColumns - 1)] = v[i];
}
return (result);
}
// . . .
public static void Test()
{
// . . .
// Example of using a homogeneous matrix to
// translate a homogeneous vector:
// Forming a homogeneous matrix that represents
// a translation:
VectorF64 homogeneousTranslationVector4x1 =
new VectorF64( 10.0, 20.0, 30.0, 1.0 );
MatrixF64 translation4x4 =
MatrixF64.FormHomogeneousTranslationMatrix
( homogeneousTranslationVector4x1 );
translation4x4.WriteLine();
// [ 1, 0, 0, 10 ]
// [ 0, 1, 0, 20 ]
// [ 0, 0, 1, 30 ]
// [ 0, 0, 0, 1 ]
// Using the homogeneous matrix to translate
// a homogeneous vector:
VectorF64 homogeneousVector =
new VectorF64( 5.0, 6.0, 7.0, 1.0 );
VectorF64 translatedVector =
translation4x4 * homogeneousVector;
translatedVector.WriteLine();
// ( 15, 26, 37, 1 )
// The relevant vector components (5,6,7)
// were translated by (10,20,30) by the
// homogeneous translation matrix.
// . . .
}
}
3.21 (d)-dimensional space : rotations
"rotation in (d)-dimensional space" : An operation (R) specified by a two-dimensional rotation plane (S), a rotation angle (t), and a rotation center point (c), that converts any point (p) to a corresponding point (q), such that:
(1) Length( q - c ) = Length( p - c ); the distance between the new point and the rotation center point is the same as the distance between the original point and the rotation center point.
Thus, the new point is on the surface of a particular (d)-dimensional sphere;
(2) Let (P) be a plane that is parallel to the rotation plane (S) and contains point (p).
Point (q) will be in plane (P).
Thus, the new point is on a particular two-dimensional plane in (d)-dimensional space;
(3) Dot( (q - c), (p - c) ) = Length( q - c ) * Length( p - c ) * Cos( t ); the dot product of the vector from the rotation center point to the original point and the vector from the rotation center point to the new point is equal to the product of the lengths of those two vectors and the cosine of the rotation angle.
Thus, the new point is on the surface of a particular (d)-dimensional cone;
Constraint (1) limits the new point to the surface of a particular (d)-dimensional sphere.
Condition (2) further constrains the new point to a particular two-dimensional plane in (d)-dimensional space.
Therefore, the new point must be on a particular two-dimensional circle in (d)-dimensional space.
Condition (3) further constrains the new point to a particular (d)-dimensional cone.
Therefore, the new point must be at one particular point (for Cos( t ) = 1 or Cos( t ) = (-1)) or at one of two possible points (for (-1) < Cos( t ) < 1).
One final constraint is required to choose one of the two possible points that result from the constraints mentioned above.
Let T(s) represents a "translation" operation, where (s) is an arbitrary vector, such that T(s) * (r) = (r + s), and Inverse(T(s)) * (r) = (r - s).
A rotation R(S,t,c), for plane S, angle t, and center point c, can be expressed in terms of translation operations and a rotation with a center point at the origin: R(S,t,c) = T(c) * R(S,t,0) * T(-c) = T(c) * R(S,t,0) * Inverse(T(c)).
Let R(S,t) represents a rotation operation with a rotation center point at the origin of (d)-dimensional space.
Thus, any arbitrary rotation can be expressed as: R(S,t,c) = T(c) * R(S,t) * T(-c) = T(c) * R(S,t) * Inverse(T(c)).
Let us consider two-dimensional rotation planes (S) that include the origin of (d)-dimensional space and are parallel to two coordinate axes.
Let (a) and (b) be integers that represent coordinate axes in (d)-dimensional space, such that 0 <= a <= (d-1), and 0 <= b <= (d-1), and (a != b).
We define R(a,b,t) to be a rotation about the origin of (d)-dimensional space, with a rotation plane parallel to coordinate axes indicated by coordinate indices (a) and (b), and with an angle (t), such that R(a,b,(pi/2)) would convert a positive coordinate value for the coordinate axis indicated by (a) to a positive coordinate value for the coordinate axis indicated by (b); example: R(a,b,(pi/2)) * e(a) = e(b), where e(k) indicates a unit basis vector corresponding to coordinate axis (k).
This definition provides the final disambiguation necessary to fully determine the new point produced by a rotation.
A rotation that moves points along axis a toward axis b according to an angle t is equivalent to a rotation that moves points along axis b toward axis a by an angle (-t); i.e., R(a,b,t) = R(b,a,(-t)).
If (d >= 2), then (d)-dimensional space has ((d*(d-1))/2) pairs of coordinate axes.
Thus, for (d >= 2), (d)-dimensional space has ((d*(d-1))/2) distinct axis-aligned rotation planes.
( Examples: d=2: { x-y plane }; d=3: { x-y plane, y-z plane, z-x plane }; d=4: { x-y plane, y-z plane, z-x plane, x-w plane, z-w plane, w-y plane }. )
The code below shows how a matrix can be used to represent the rotation R(a,b,t).
An example is shown that uses a rotation matrix in a homogeneous manner to rotate a homogeneous vector to a new homogeneous vector.
public class MatrixF64
{
// . . .
public static MatrixF64 Rab( int rows, int a, int b, double angle )
{
if (rows <= 0)
{
// Invalid row count specified.
return (MatrixF64.Zero( 1, 1 ));
}
MatrixF64 result = MatrixF64.Identity( rows );
// If (a==b) (rotate coordinate (a) to (a)), then
// there is no rotation to do.
if ((a != b) && (a >= 0) && (a < rows) && (b >= 0) && (b < rows))
{
double sine = Math.Sin( angle );
double cosine = Math.Cos( angle );
// Override m[a,a] = cosine; m[a,b] = -sine;
// m[b,a] = sine; m[b,b] = cosine;
result[a, a] = cosine;
result[b, b] = cosine;
result[a, b] = (-(sine));
result[b, a] = sine;
}
return (result);
}
public static MatrixF64 Rabk( int rows, int a, int b, int k )
{
double angle = (Math.PI / 2.0) * (double)k;
return (MatrixF64.Rab( rows, a, b, angle ));
}
// . . .
public static void Test()
{
// . . .
// Example of using a homogeneous matrix
// to rotate a homogeneous vector:
// Forming a homogeneous rotation matrix that
// rotates a positive component along coordinate
// axis 0 to a positive component along
// cooridnate axis 1, in 3-dimensional space.
// (The homogeneous matrix must be (3+1)=4
// rows by (3+1)=4 columns.)
MatrixF64 rotation4x4 =
MatrixF64.Rab( 4, 0, 1, (0.5 * Math.PI) );
rotation4x4.WriteLine(3);
// [ 6.12e-17, -1, 0, 0 ]
// [ 1, 6.12e-17, 0, 0 ]
// [ 0, 0, 1, 0 ]
// [ 0, 0, 0, 1 ]
// (The entries '6.12e-17' are negligible.)
// The upper-left 3x3 part of the matrix represents
// a counterclockwise rotation of a quarter-turn
// such that a positive component along axis 0
// would be rotated to a positive component along
// axis 1, in 3-dimensional space.
// Using the homogeneous rotation matrix to
// rotate a homogeneous vector:
VectorF64 homogeneousVector2 =
new VectorF64( 5.0, 6.0, 7.0, 1.0 );
VectorF64 rotatedVector =
rotation4x4 * homogeneousVector2;
rotatedVector.WriteLine();
// ( -6, 5, 7, 1 )
// The relevant vector components (5,6,7) were
// rotated to (-6,5,7), which is consistent with
// a counterclockwise quarter-turn rotation
// with a rotation plane parallel to axes 0 and 1.
// . . .
}
}
3.22 Proposed definition of "counterclockwise"
In two-dimensional space, the quarter-turn rotation that can convert the vector e(0) to e(1) is considered counterclockwise.
The rotation is: R(0,1,(pi/2)) * e(0) = e(1).
Therefore, the rotation R(0,1,t) is counterclockwise.
In three-dimensional space, the quarter-turn rotation that can convert the vector e(0) to e(1) is considered counterclockwise.
The quarter-turn rotation that can convert the vector e(1) to e(2) is also considered counterclockwise.
The quarter-turn rotation that can convert the vector e(2) to e(0) is also considered counterclockwise.
Therefore, the rotations { R(0,1,t), R(1,2,t), R(2,0,t) } are counterclockwise.
The Boolean formula below fits the definition of counterclockwise for two-dimensional and three-dimensional space, and extrapolates the pattern for the general case of the rotation R(a,b,t):
bool counterClockwise =
(
( (a < b) && (((a + b) % 2) == 1) )
|| ( (a > b) && (((a + b) % 2) == 0) )
);
This rule would define the rotations below as "counterclockwise":
R(0,1,t),
R(2,0,t), R(1,2,t),
R(0,3,t), R(3,1,t), R(2,3,t),
R(4,0,t), R(1,4,t), R(4,2,t), R(3,4,t),
R(0,5,t), R(5,1,t), R(2,5,t), R(5,3,t), R(4,5,t),
...
3.23 (d)-dimensional space : clockwise rotations
Using the proposed definition of "counterclockwise" in the previous section, the rotations below are "clockwise":
R(1,0,t),
R(0,2,t), R(2,1,t),
R(3,0,t), R(1,3,t), R(3,2,t),
R(0,4,t), R(4,1,t), R(2,4,t), R(4,3,t),
R(5,0,t), R(1,5,t), R(5,2,t), R(3,5,t), R(5,4,t),
...
For example, in four-dimensional space, the six rotation matrices below represent clockwise rotations:
R(1,0,t) = Cos(t), Sin(t), 0, 0,
-Sin(t), Cos(t), 0, 0,
0, 0, 1, 0,
0, 0, 0, 1
R(0,2,t) = Cos(t), 0, -Sin(t), 0,
0, 1, 0, 0,
Sin(t), 0, Cos(t), 0,
0, 0, 0, 1
R(2,1,t) = 1, 0, 0, 0,
0, Cos(t), Sin(t), 0,
0, -Sin(t), Cos(t), 0,
0, 0, 0, 1
R(3,0,t) = Cos(t), 0, 0, Sin(t),
0, 1, 0, 0,
0, 0, 1, 0,
-Sin(t), 0, 0, Cos(t)
R(1,3,t) = 1, 0, 0, 0,
0, Cos(t), 0, -Sin(t),
0, 0, 1, 0,
0, Sin(t), 0, Cos(t)
R(3,2,t) = 1, 0, 0, 0,
0, 1, 0, 0,
0, 0, Cos(t), Sin(t),
0, 0, -Sin(t), Cos(t)
3.24 (d)-dimensional space : total number of axis-aligned orientations
The number of distinct axis-aligned orientations in (d)-dimensional space is:
distinctOrientations = Factorial( d ) * IntPow( 2, (d-1) );
An argument supporting this formula appears in the paragraphs below.
Experimental evidence verifying the formula is also presented below.
Let { g(0), g(1), ..., g(d-1) } be a set of (d) mutually-perpendicular unit vectors attached to an object in (d)-dimensional space.
These unit vectors indicate the orientation of the object relative to the coordinate axes of the (d)-dimensional space.
Let the "default orientation" of the object be such that g(i) = e(i), for all 0 <= i <= (d-1); each orientation vector of the object is associated with a distinct unit basis vector of a coordinate system in the (d)-dimensional space.
The first orientation vector of the object can be parallel with any of the (d) axes of the (d)-dimensional space, and can be in either of two directions relative to that axis.
Thus, the first orientation vector has (2 * d) options.
The second orientation vector of the object can be parallel to any of the (d-1) remaining axes, with two possible directions relative to that axis.
Thus, the second orientation vector has (2 * (d-1)) options.
In the same way, the next orientation vector has (2 * (d-2)) options.
The final orientation vector is forced to the final remaining axis, and is forced to one of the two possible directions along that axis (to conserve the parity of the set of orientation vectors; the determinant of an ordered set of the orientation vectors must be constant).
Thus, the final orientation vector only has (1) option.
The total number of combinations is given by the product: ((2 * d) * (2 * (d-1)) * ... * (2 * 3) * (2 * 2) * (1)) = (2^(d-1)) * (d*(d-1)*(d-2)*...*2*1) = (2^(d-1)) * Factorial(d).
The paragraphs below present reassuring experimental evidence supporting the formula for the total number of axis-aligned orientations in which an object can be.
The orientation of an object can be represented by a matrix (M) where the vectors { g(0), g(1), ..., g(d-1) } are the columns of the matrix:
M = g(0)[ 0 ], g(1)[ 0 ], ..., g(d-2)[ 0 ], g(d-1)[ 0 ],
g(0)[ 1 ], g(1)[ 1 ], ..., g(d-2)[ 1 ], g(d-1)[ 1 ],
g(0)[ 2 ], g(1)[ 2 ], ..., g(d-2)[ 2 ], g(d-1)[ 2 ],
g(0)[ 3 ], g(1)[ 3 ], ..., g(d-2)[ 3 ], g(d-1)[ 3 ],
. . . . .
. . . . .
g(0)[d-4], g(1)[d-4], ..., g(d-2)[d-4], g(d-1)[d-4],
g(0)[d-3], g(1)[d-3], ..., g(d-2)[d-3], g(d-1)[d-3],
g(0)[d-2], g(1)[d-2], ..., g(d-2)[d-2], g(d-1)[d-2],
g(0)[d-1], g(1)[d-1], ..., g(d-2)[d-1], g(d-1)[d-1]
The matrix (M) corresponding to the "default orientation" is the "identity matrix":
M = 1, 0, ..., 0, 0,
0, 1, ..., 0, 0,
0, 0, ..., 0, 0,
0, 0, ..., 0, 0,
0, 0, ..., 0, 0,
. . . . .
. . . . .
0, 0, ..., 0, 0,
0, 0, ..., 0, 0,
0, 0, ..., 1, 0,
0, 0, ..., 0, 1
A rotation matrix, R(a,b,(pi/2)), corresponding to a quarter-turn rotation such that e(a) is rotated to e(b), can convert the orientation matrix of the object to a new orientation matrix; R(a,b,(pi/2)) * M = Mnew.
For (d) dimensions, we consider all ((d*(d-1))/2) combinations of coordinate axes (a) and (b), and the corresponding quarter-turn rotations R(a,b,(pi/2)).
For example, in 4-dimensional space, we can consider the ((4*(4-1))/2) = 6 rotations in the set { R(1,0,(pi/2)), R(0,2,(pi/2)), R(2,1,(pi/2)), R(3,0,(pi/2)), R(1,4,(pi/2)), R(3,2,(pi/2)) }.
( Only clockwise rotations need be considered, because combining clockwise rotations can create counterclockwise rotations. )
The total number of axis-aligned orientations that an object can be in can be determined by starting with an initial orientation (represented by the identity matrix), then applying each possible quarter-turn matrix to the initial matrix to find new orientation matrices, and then applying each possible quarter-turn matrix to all found orientation matrices to find more distinct orientation matrices.
Continue until no new matrices can be found, such that applying all possible quarter-turn matrices to all found matrices yields no new matrices.
The number of distinct matrices found by this process is the number of distinct axis-aligned orientations an object can be in.
An exhaustive search produces the results below:
d = 2: 4 axis-aligned orientations
d = 3: 24 axis-aligned orientations
d = 4: 192 axis-aligned orientations
d = 5: 1920 axis-aligned orientations
d = 6: 23040 axis-aligned orientations
d = 7: 322560 axis-aligned orientations
d = 8: 5160960 axis-aligned orientations
. . . . . . . . . .
These experimental results are consistent with the numbers produced by the general formula.
Only two of the three quarter-turn rotations in the set { R(1,0,(pi/2)), R(0,2,(pi/2)), R(2,1,(pi/2)) } are required to rotate an object in to all 24 possible axis-aligned orientations in three-dimensional space.
We can, for example, do without R(1,0,(pi/2)), because:
R(1,0,(pi/2)) = (R(2,1,(pi/2)))^3 * R(0,2,(pi/2) * R(2,1,(pi/2));
where,
(R(2,1,(pi/2)))^3 = R(2,1,(pi/2)) * R(2,1,(pi/2)) * R(2,1,(pi/2));
For example, if { R(1,0,(pi/2)), R(0,2,(pi/2)), R(2,1,(pi/2)) } are clockwise quarter-turn rotations about the { Z, Y, X } axes, then it is easy to verify, by holding and turning a physical object, that doing clockwise quarter-turns, first about the X axis [see note below], then about the Y axis, and then three quarter-turns about the X axis, will result in an orientation that is equivalent to doing a clockwise quarter-turn rotation about the Z axis.
[ Note: The right-most operator in any mathematical product of operators is the operator with the earliest effect on an object to which the overall product is applied.
For example, given operators (A), (B), (C), the product (C * B * A) means that (A) has the earliest effect.
This can be seen by considering a particular object (P) and noting that the associative attribute of the operators under consideration implies (C * B * A) * P = (C * B) * (A * P) = C * (B * (A * P)).
Thus, the mathematical product "(R(2,1,(pi/2)))^3 * R(0,2,(pi/2) * R(2,1,(pi/2))" implies that "R(2,1,(pi/2))" has the earliest effect, and "R(0,2,(pi/2)" has the next effect, and "(R(2,1,(pi/2)))^3" acts last. ]
Therefore, it would be sufficient, for example, to use only the rotations { R(0,2,(pi/2)), R(2,1,(pi/2)) } to access all 24 axis-aligned orientations in three-dimensional space.
Likewise, only three of the six quarter-turn rotations in the set { R(1,0,(pi/2)), R(0,2,(pi/2)), R(2,1,(pi/2)), R(3,0,(pi/2)), R(1,4,(pi/2)), R(3,2,(pi/2)) } are required to rotate and object in to all 192 possible axis-aligned orientations in four-dimensional space.
We can, for example, do without R(1,0,(pi/2)), R(3,2,(pi/2)), and R(1,3,(pi/2)), because:
R(1,0,(pi/2)) = (R(2,1,(pi/2)))^3 * R(0,2,(pi/2)) * R(2,1,(pi/2));
R(3,2,(pi/2)) = (R(3,0,(pi/2)))^3 * R(0,2,(pi/2)) * R(3,0,(pi/2));
R(1,3,(pi/2)) = (R(3,0,(pi/2)))^3 * R(1,0,(pi/2)) * R(3,0,(pi/2));
Therefore, it would be sufficient, for example, to use only the rotations { R(0,2,(pi/2)), R(2,1,(pi/2)), R(3,0,(pi/2)) } to access all 192 axis-aligned orientations in four-dimensional space.
3.25 3-dimensional space and 4-dimensional space : axis-aligned orientations
In 3-dimensional space there are 24 axis-aligned orientations.
Each orientation is linked to 9 other orientations by the quarter-turn, half-turn, and inverse quarter-turn, rotations below:
Z = R(1,0,(pi/2)), Z2 = R(1,0,pi), Zi = R(1,0,(-(pi/2))),
Y = R(0,2,(pi/2)), Y2 = R(0,2,pi), Yi = R(0,2,(-(pi/2))),
X = R(2,1,(pi/2)), X2 = R(2,1,pi), Xi = R(2,1,(-(pi/2)))
The symbols { Z, Z2, Zi, Y, Y2, Yi, X, X2, Xi } are abbreviations for the rotations.
The diagram below shows how each axis-aligned orientation is directly linked to 9 other orientations by the 9 rotations { Z, Z2, Zi, Y, Y2, Yi, X, X2, Xi }.
The table below shows how each of the 24 axis-aligned orientations is directly linked to 9 other axis-aligned orientations by the 9 rotations { Z, Z2, Zi, Y, Y2, Yi, X, X2, Xi }.
Each orientation is indicated by three orientation vectors, { G0, G1, G2 }, where each axis-aligned orientation vector can point in one of six possible directions, { -X, X, -Y, Y, -Z, Z }, relative to reference coordinate system axes.
The table entries show the orientation (1..24) that results from applying the rotation operation associated with a column of the table to the orientation associated with a row of the table.
-------------------------------------------------------------
# : G0 G1 G2 : X X2 Xi Y Y2 Yi Z Z2 Zi
-------------------------------------------------------------
1 : X Y Z : 5 2 6 14 3 16 24 4 21
2 : X -Y -Z : 6 1 5 13 4 15 23 3 22
3 : -X Y -Z : 7 4 8 16 1 14 22 2 23
4 : -X -Y Z : 8 3 7 15 2 13 21 1 24
-------------------------------------------------------------
5 : X -Z Y : 2 6 1 9 8 12 20 7 18
6 : X Z -Y : 1 5 2 10 7 11 19 8 17
7 : -X -Z -Y : 4 8 3 11 6 10 18 5 20
8 : -X Z Y : 3 7 4 12 5 9 17 6 19
-------------------------------------------------------------
9 : Z X Y : 22 11 24 8 12 5 13 10 14
10 : Z -X -Y : 21 12 23 7 11 6 14 9 13
11 : -Z X -Y : 24 9 22 6 10 7 15 12 16
12 : -Z -X Y : 23 10 21 5 9 8 16 11 15
-------------------------------------------------------------
13 : Z -Y X : 17 16 20 4 15 2 10 14 9
14 : Z Y -X : 18 15 19 3 16 1 9 13 10
15 : -Z -Y -X : 19 14 18 2 13 4 12 16 11
16 : -Z Y X : 20 13 17 1 14 3 11 15 12
-------------------------------------------------------------
17 : Y Z X : 16 20 13 21 18 22 6 19 8
18 : Y -Z -X : 15 19 14 22 17 21 5 20 7
19 : -Y Z -X : 14 18 15 23 20 24 8 17 6
20 : -Y -Z X : 13 17 16 24 19 23 7 18 5
-------------------------------------------------------------
21 : Y -X Z : 12 23 10 18 22 17 1 24 4
22 : Y X -Z : 11 24 9 17 21 18 2 23 3
23 : -Y -X -Z : 10 21 12 20 24 19 3 22 2
24 : -Y X Z : 9 22 11 19 23 20 4 21 1
-------------------------------------------------------------
If we permit all 9 rotations from the set { Z, Z2, Zi, Y, Y2, Yi, X, X2, Xi }, then the resulting histogram of the number of orientations requiring a minimum of a certain number of rotations to reach is: { 0:1, 1:9, 2:14 };
a minimum of 0 rotations is required to reach 1 orientation (the given orientation itself), and a minimum of 1 rotation is required to reach 9 orientations (where each orientation is produced by rotating the given orientation by one of the 9 permitted rotations), and a minimum of 2 rotations is required to reach 14 orientations.
The histogram indicates that permitting the 9 rotations makes it possible to reach all 1 + 9 + 14 = 24 orientations.
The histogram indicates that the number of rotations, from the set of permitted rotations, required to convert any given orientation to any other given orientation, will be two or fewer rotations.
If the set of permitted rotations is arbitrarily reduced, then the number of links between orientations will be reduced, and the minimum number of rotations required to convert a given orientation to another given orientation might be larger than when more rotations are permitted.
If we permit all 9 rotations from the set { Z, Z2, Zi, Y, Y2, Yi, X, X2, Xi }, then the resulting histogram of the number of orientations requiring a minimum of a certain number of rotations to reach is: { 0:1, 1:9, 2:14 }.
If we only permit the 3 rotations from the set { Z, Y, X } (only the clockwise quarter-turn rotations), then, for any given orientation, the histogram becomes: { 0:1, 1:6, 2:11, 3:6 }, which indicates that all 1 + 6 + 11 + 6 = 24 orientations are reachable from the given orientation, and that three or fewer rotations are required to convert any given orientation to any other given orientation.
If we only permit the 2 rotations from the set { Z, Y }, then, for any given orientation, the histogram becomes: { 0:1, 1:4, 2:10, 3:8, 4:1 }, which indicates that all 1 + 4 + 10 + 8 + 1 = 24 orientations are reachable from the given orientation, and that four or fewer rotations are required to convert any given orientation to any other given orientation.
( X = (Z * Y * Zi) = (Z * Y * Z*Z*Z).
Therefore, { Z, Y } is sufficient to reach all orientations that can be reached with { Z, Y, X } (all 24 orientations).)
If we only permit the 1 rotation from the set { Z }, then, for any given orientation, the histogram becomes: { 0:1, 1:2, 2:1 }, which indicates that only 1 + 2 + 1 = 4 orientations are reachable from the given orientation; and, therefore, the set of 24 orientations is broken in to (24/4) = 6 distinct subsets of orientations that contain mutually-accessible orientations.
In three-dimensional space, a minimum of two rotations (associated with distinct pairs of axes, and satisfying other conditions) are required to make it possible to convert any orientation to all other orientations.
In 4-dimensional space there are 192 axis-aligned orientations.
Each orientation is linked to 18 other orientations by the quarter-turn, half-turn, and inverse quarter-turn, rotations below:
Z = R(1,0,(pi/2)), Z2 = R(1,0,pi), Zi = R(1,0,(-(pi/2))),
Y = R(0,2,(pi/2)), Y2 = R(0,2,pi), Yi = R(0,2,(-(pi/2))),
X = R(2,1,(pi/2)), X2 = R(2,1,pi), Xi = R(2,1,(-(pi/2))),
WX = R(3,0,(pi/2)), WX2 = R(3,0,pi), WXi = R(3,0,(-(pi/2))),
YW = R(1,3,(pi/2)), YW2 = R(1,3,pi), YWi = R(1,3,(-(pi/2))),
WZ = R(3,2,(pi/2)), WZ2 = R(3,2,pi), WZi = R(3,2,(-(pi/2)))
The symbols { Z, Z2, Zi, Y, Y2, Yi, X, X2, Xi, WX, WX2, WXi, YW, YW2, YWi, WZ, WZ2, WZi } are abbreviations for the rotations.
The diagram below shows how each axis-aligned orientation is directly linked to 18 other orientations by the 18 rotations { Z, Z2, Zi, Y, Y2, Yi, X, X2, Xi, WX, WX2, WXi, YW, YW2, YWi, WZ, WZ2, WZi }.
Each of the 192 axis-aligned orientations is directly linked to 18 other axis-aligned orientations by the 18 rotations { Z, Z2, Zi, Y, Y2, Yi, X, X2, Xi, WX, WX2, WXi, YW, YW2, YWi, WZ, WZ2, WZi }.
Each orientation can be indicated by four orientation vectors, { G0, G1, G2, G3 }, where each axis-aligned orientation vector can point in one of eight possible directions, { -X, X, -Y, Y, -Z, Z, -W, W }, relative to reference coordinate system axes.
If we permit all 18 rotations from the set { Z, Z2, Zi, Y, Y2, Yi, X, X2, Xi, WX, WX2, WXi, YW, YW2, YWi, WZ, WZ2, WZi }, then the resulting histogram of the number of orientations requiring a minimum of a certain number of rotations to reach is:
{ 0:1, 1:18, 2:81, 3:92 };
a minimum of 0 rotations is required to reach 1 orientation (the given orientation itself),
and a minimum of 1 rotation is required to reach 18 orientations (where each orientation is produced by rotating the given orientation by one of the 18 permitted rotations),
and a minimum of 2 rotations is required to reach 81 orientations,
and a minimum of 3 rotations is required to reach 92 orientations.
The histogram indicates that permitting the 18 rotations makes it possible to reach all 1 + 18 + 81 + 92 = 192 orientations.
The histogram indicates that the number of rotations, from the set of permitted rotations, required to convert any given orientation to any other given orientation, will be three or fewer rotations.
If the set of permitted rotations is arbitrarily reduced, then the number of links between orientations will be reduced, and the minimum number of rotations required to convert a given orientation to another given orientation might be larger than when more rotations are permitted.
If we permit all 18 rotations from the set { Z, Z2, Zi, Y, Y2, Yi, X, X2, Xi, WX, WX2, WXi, YW, YW2, YWi, WZ, WZ2, WZi }, then the resulting histogram of the number of orientations requiring a minimum of a certain number of rotations to reach is: { 0:1, 1:18, 2:81, 3:92 }.
If we only permit the 6 rotations from the set { Z, Y, X, WX, YW, WZ } (only the clockwise quarter-turn rotations), then, for any given orientation, the histogram becomes: { 0:1, 1:12, 2:50, 3:84, 4:45 }, which indicates that all 1 + 12 + 50 + 84 + 45 = 192 orientations are reachable from the given orientation, and that four or fewer rotations are required to convert any given orientation to any other given orientation.
If we only permit the 3 rotations from the set { Z, Y, WX }, then, for any given orientation, the histogram becomes:
{ 0:1, 1:6, 2:27, 3:72, 4:67, 5:18, 6:1 }, which indicates that all 1 + 6 + 27 + 72 + 67 + 18 + 1 = 192 orientations are reachable from the given orientation, and that six or fewer rotations are required to convert any given orientation to any other given orientation.
If we only permit rotations from the set { Z, X, WX }, or only from the set { Y, X, WX }, then the histogram becomes: { 0:1, 1:6, 2:23, 3:52, 4:63, 5:38, 6:9 }, which indicates that all 192 orientations are accessible, but that more rotations are required to reach some orientations than would be if we instead only permitted rotations from the set { Z, Y, WX }.
If we only permit the 2 rotations from the set { Z, Y }, then, for any given orientation, the histogram becomes: { 0:1, 1:4, 2:10, 3:8, 4:1 }, which indicates that only 1 + 4 + 10 + 8 + 1 = 24 orientations are reachable from the given orientation;
and, therefore, the set of 192 orientations is broken in to (192/24) = 8 distinct subsets of orientations that contain mutually-accessible orientations.
In four-dimensional space, a minimum of three rotations (associated with distinct pairs of axes, and satisfying other conditions) are required to make it possible to convert any orientation to all other orientations.
