This file provides linear algebra, with support for matrix construction, solutions to linear systems of equations, QR and Cholesky factorizations, and matrix inverse.

To use, add the following lines to the beginning of your file:

include <BOSL2/std.scad>

1. Section: Matrices

2. Section: Matrix testing and display

   • is_matrix() – Check if input is a numeric matrix, optionally of specified size
   • is_matrix_symmetric() – Checks if matrix is symmetric
   • is_rotation() – Check if a transformation matrix represents a rotation.
   • echo_matrix() – Print a matrix neatly to the console.

3. Section: Matrix indexing

   • column() – Extract a column from a matrix.
   • submatrix() – Extract a submatrix from a matrix

4. Section: Matrix construction and modification

   • ident() – Return identity matrix.
   • diagonal_matrix() – Make a diagonal matrix.
   • transpose() – Transpose a matrix
   • outer_product() – Compute the outer product of two vectors.
   • submatrix_set() – Takes a matrix as input and change values in a submatrix.
   • hstack() – Make a new matrix by stacking matrices horizontally.
   • block_matrix() – Make a new matrix from a block of matrices.

5. Section: Solving Linear Equations and Matrix Factorizations

   • linear_solve() – Solve Ax=b or, for overdetermined case, solve the least square problem.
   • linear_solve3() – Fast solution to Ax=b where A is 3x3.
   • matrix_inverse() – General matrix inverse.
   • rot_inverse() – Invert 2d or 3d rotation transformations.
   • null_space() – Return basis for the null space of A.
   • qr_factor() – Compute QR factorization of a matrix.
   • back_substitute() – Solve an upper triangular system, Rx=b.
   • cholesky() – Compute the Cholesky factorization of a matrix.

6. Section: Matrix Properties: Determinants, Norm, Trace

   • det2() – Compute determinant of 2x2 matrix.
   • det3() – Compute determinant of 3x3 matrix.
   • det4() – Compute determinant of 4x4 matrix.
   • determinant() – compute determinant of an arbitrary square matrix.
   • norm_fro() – Compute Frobenius norm of a matrix
   • matrix_trace() – Compute the trace of a square matrix.

The matrix, a rectangular array of numbers which represents a linear transformation, is the fundamental object in linear algebra. In OpenSCAD a matrix is a list of lists of numbers with a rectangular structure. Because OpenSCAD treats all data the same, most of the functions that index matrices or construct them will work on matrices (lists of lists) whose elements are not numbers but may be arbitrary data: strings, booleans, or even other lists. It may even be acceptable in some cases if the structure is non-rectangular. Of course, linear algebra computations and solutions require true matrices with rectangular structure, where all the entries are finite numbers.

Matrices in OpenSCAD are lists of row vectors. However, a potential source of confusion is that OpenSCAD treats vectors as either column vectors or row vectors as demanded by context. Thus both v*M and M*v are valid if M is square and v has the right length. If you want to multiply M on the left by v and w you can do this with [v,w]*M but if you want to multiply on the right side with v and w as column vectors, you now need to use transpose() because OpenSCAD doesn't adjust matrices contextually: A=M*transpose([v,w]). The solutions are now columns of A and you must extract them with column() or take the transpose of A.

Synopsis: Check if input is a numeric matrix, optionally of specified size

Topics: Matrices

See Also: is_matrix_symmetric(), is_rotation()

Usage:

• test = is_matrix(A, [m], [n], [square])

Description:

Returns true if A is a numeric matrix of height m and width n with finite entries. If m or n are omitted or set to undef then true is returned for any positive dimension.

Arguments:

Synopsis: Checks if matrix is symmetric

Topics: Matrices

See Also: is_matrix(), is_rotation()

Usage:

• b = is_matrix_symmetric(A, [eps])

Description:

Returns true if the input matrix is symmetric, meaning it approximately equals its transpose. The matrix can have arbitrary entries.

Arguments:

Synopsis: Check if a transformation matrix represents a rotation.

Topics: Affine, Matrices, Transforms

See Also: is_matrix(), is_matrix_symmetric()

Usage:

• b = is_rotation(A, [dim], [centered])

Description:

Returns true if the input matrix is a square affine matrix that is a rotation around any point, or around the origin if centered is true. The matrix must be 3x3 (representing a 2d transformation) or 4x4 (representing a 3d transformation). You can set dim to 2 to require a 2d transform (3x3 matrix) or to 3 to require a 3d transform (4x4 matrix).

Arguments:

Synopsis: Print a matrix neatly to the console.

Topics: Matrices

See Also: is_matrix(), is_matrix_symmetric(), is_rotation()

Usage:

• echo_matrix(M, [description], [sig], [sep], [eps]);
• dummy = echo_matrix(M, [description], [sig], [sep], [eps]),

Description:

Display a numerical matrix in a readable columnar format with sig significant digits. Values smaller than eps display as zero. If you give a description it is displayed at the top. You can change the space between columns by setting sep to a number of spaces, which will use wide figure spaces the same width as digits, or you can set it to any string to separate the columns. Values that are NaN or INF will display as "nan" and "inf". Values which are otherwise non-numerica display as two dashes. Note that this includes lists, so a 3D array will display as a list of dashes.

Arguments:

Synopsis: Extract a column from a matrix.

Topics: Matrices, List Handling, Arrays

See Also: select(), slice()

Usage:

• list = column(M, i);

Description:

Extracts entry i from each list in M, or equivalently column i from the matrix M, and returns it as a vector. This function will return undef at all entry positions indexed by i not found in M.

Arguments:

Example 1:

include <BOSL2/std.scad>
M = [[1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]];
a = column(M,2);      // Returns [3, 7, 11, 15]
b = column(M,0);      // Returns [1, 5, 9, 13]
N = [ [1,2], [3], [4,5], [6,7,8] ];
c = column(N,1);      // Returns [1,undef,5,7]
data = [[1,[3,4]], [3, [9,3]], [4, [3,1]]];   // Matrix with non-numeric entries
d = column(data,0);   // Returns [1,3,4]
e = column(data,1);   // Returns [[3,4],[9,3],[3,1]]

Synopsis: Extract a submatrix from a matrix

Topics: Matrices, Arrays

See Also: column(), block_matrix(), submatrix_set()

Usage:

• mat = submatrix(M, idx1, idx2);

Description:

The input must be a list of lists (a matrix or 2d array). Returns a submatrix by selecting the rows listed in idx1 and columns listed in idx2.

Arguments:

Example 1:

include <BOSL2/std.scad>
M = [[ 1, 2, 3, 4, 5],
     [ 6, 7, 8, 9,10],
     [11,12,13,14,15],
     [16,17,18,19,20],
     [21,22,23,24,25]];
submatrix(M,[1:2],[3:4]);  // Returns [[9, 10], [14, 15]]
submatrix(M,[1], [3,4]));  // Returns [[9,10]]
submatrix(M,1, [3,4]));  // Returns [[9,10]]
submatrix(M,1,3));  // Returns [[9]]
submatrix(M, [3,4],1); // Returns  [[17],[22]]);
submatrix(M, [1,3],[2,4]); // Returns [[8,10],[18,20]]);
A = [[true,    17, "test"],
     [[4,2],   91, false],
     [6,    [3,4], undef]];
submatrix(A,[0,2],[1,2]);   // Returns [[17, "test"], [[3, 4], undef]]

Synopsis: Return identity matrix.

Topics: Affine, Matrices, Transforms

See Also: IDENT, submatrix(), column()

Usage:

• mat = ident(n);

Description:

Create an n by n square identity matrix.

Arguments:

Example 1:

include <BOSL2/std.scad>
mat = ident(3);
// Returns:
//   [
//     [1, 0, 0],
//     [0, 1, 0],
//     [0, 0, 1]
//   ]

Example 2:

include <BOSL2/std.scad>
mat = ident(4);
// Returns:
//   [
//     [1, 0, 0, 0],
//     [0, 1, 0, 0],
//     [0, 0, 1, 0],
//     [0, 0, 0, 1]
//   ]

Synopsis: Make a diagonal matrix.

Topics: Affine, Matrices

See Also: column(), submatrix()

Usage:

• mat = diagonal_matrix(diag, [offdiag]);

Description:

Creates a square matrix with the items in the list diag on its diagonal. The off diagonal entries are set to offdiag, which is zero by default.

Arguments:

Synopsis: Transpose a matrix

Topics: Linear Algebra, Matrices

See Also: submatrix(), block_matrix(), hstack(), flatten()

Usage:

• M = transpose(M, [reverse]);

Description:

Returns the transpose of the given input matrix. The input can be a matrix with arbitrary entries or a numerical vector. If you give a vector then transpose returns it unchanged. When reverse=true, the transpose is done across to the secondary diagonal. (See example below.) By default, reverse=false.

Example 1:

include <BOSL2/std.scad>
M = [
    [1, 2, 3],
    [4, 5, 6],
    [7, 8, 9]
];
t = transpose(M);
// Returns:
// [
//     [1, 4, 7],
//     [2, 5, 8],
//     [3, 6, 9]
// ]

Example 2:

include <BOSL2/std.scad>
M = [
    [1, 2, 3],
    [4, 5, 6]
];
t = transpose(M);
// Returns:
// [
//     [1, 4],
//     [2, 5],
//     [3, 6],
// ]

Example 3:

include <BOSL2/std.scad>
M = [
    [1, 2, 3],
    [4, 5, 6],
    [7, 8, 9]
];
t = transpose(M, reverse=true);
// Returns:
// [
//  [9, 6, 3],
//  [8, 5, 2],
//  [7, 4, 1]
// ]

Example 4: Transpose on a list of numbers returns the list unchanged

include <BOSL2/std.scad>
transpose([3,4,5]);  // Returns: [3,4,5]

Example 5: Transpose on non-numeric input

include <BOSL2/std.scad>
arr = [
    [  "a",  "b", "c"],
    [  "d",  "e", "f"],
    [[1,2],[3,4],[5,6]]
];
t = transpose(arr);
// Returns:
// [
//     ["a", "d", [1,2]],
//     ["b", "e", [3,4]],
//     ["c", "f", [5,6]],
// ]

Synopsis: Compute the outer product of two vectors.

Topics: Linear Algebra, Matrices

See Also: submatrix(), determinant()

Usage:

• x = outer_product(u,v);

Usage:

• M = outer_product(u,v);

Description:

Compute the outer product of two vectors, which is a matrix.

Synopsis: Takes a matrix as input and change values in a submatrix.

Topics: Matrices, Arrays

See Also: column(), submatrix()

Usage:

• mat = submatrix_set(M, A, [m], [n]);

Description:

Sets a submatrix of M equal to the matrix A. By default the top left corner of M is set to A, but you can specify offset coordinates m and n. If A (as adjusted by m and n) extends beyond the bounds of M then the extra entries are ignored. You can pass in A=[[]], a null matrix, and M will be returned unchanged. This function works on arbitrary lists of lists and the input M need not be rectangular in shape.

Arguments:

Synopsis: Make a new matrix by stacking matrices horizontally.

Topics: Matrices, Arrays

See Also: column(), submatrix(), block_matrix()

Usage:

• A = hstack(M1, M2)
• A = hstack(M1, M2, M3)
• A = hstack([M1, M2, M3, ...])

Description:

Constructs a matrix by horizontally "stacking" together compatible matrices or vectors. Vectors are treated as columsn in the stack. This command is the inverse of column. Note: strings given in vectors are broken apart into lists of characters. Strings given in matrices are preserved as strings. If you need to combine vectors of strings use list_to_matrix() as shown below to convert the vector into a column matrix. Also note that vertical stacking can be done directly with concat.

Arguments:

Example 1:

include <BOSL2/std.scad>
M = ident(3);
v1 = [2,3,4];
v2 = [5,6,7];
v3 = [8,9,10];
a = hstack(v1,v2);     // Returns [[2, 5], [3, 6], [4, 7]]
b = hstack(v1,v2,v3);  // Returns [[2, 5,  8],
                       //          [3, 6,  9],
                       //          [4, 7, 10]]
c = hstack([M,v1,M]);  // Returns [[1, 0, 0, 2, 1, 0, 0],
                       //          [0, 1, 0, 3, 0, 1, 0],
                       //          [0, 0, 1, 4, 0, 0, 1]]
d = hstack(column(M,0), submatrix(M,idx(M),[1 2]));  // Returns M
strvec = ["one","two"];
strmat = [["three","four"], ["five","six"]];
e = hstack(strvec,strvec); // Returns [["o", "n", "e", "o", "n", "e"],
                           //          ["t", "w", "o", "t", "w", "o"]]
f = hstack(list_to_matrix(strvec,1), list_to_matrix(strvec,1));
                           // Returns [["one", "one"],
                           //          ["two", "two"]]
g = hstack(strmat,strmat); //  Returns: [["three", "four", "three", "four"],
                           //            [ "five",  "six",  "five",  "six"]]

Synopsis: Make a new matrix from a block of matrices.

Topics: Matrices, Arrays

See Also: column(), submatrix()

Usage:

• bmat = block_matrix([[M11, M12,...],[M21, M22,...], ... ]);

Description:

Create a block matrix by supplying a matrix of matrices, which will be combined into one unified matrix. Every matrix in one row must have the same height, and the combined width of the matrices in each row must be equal. Strings will stay strings.

Example 1:

include <BOSL2/std.scad>
A = [[1,2],
     [3,4]];
B = ident(2);
C = block_matrix([[A,B],[B,A],[A,B]]);
    // Returns:
    //        [[1, 2, 1, 0],
    //         [3, 4, 0, 1],
    //         [1, 0, 1, 2],
    //         [0, 1, 3, 4],
    //         [1, 2, 1, 0],
    //         [3, 4, 0, 1]]);
D = block_matrix([[A,B], ident(4)]);
    // Returns:
    //        [[1, 2, 1, 0],
    //         [3, 4, 0, 1],
    //         [1, 0, 0, 0],
    //         [0, 1, 0, 0],
    //         [0, 0, 1, 0],
    //         [0, 0, 0, 1]]);
E = [["one", "two"], [3,4]];
F = block_matrix([[E,E]]);
    // Returns:
    //        [["one", "two", "one", "two"],
    //         [    3,     4,     3,     4]]

Synopsis: Solve Ax=b or, for overdetermined case, solve the least square problem.

Topics: Matrices, Linear Algebra

See Also: linear_solve3(), matrix_inverse(), rot_inverse(), back_substitute(), cholesky()

Usage:

• solv = linear_solve(A,b,[pivot])

Description:

Solves the linear system Ax=b. If A is square and non-singular the unique solution is returned. If A is overdetermined the least squares solution is returned. If A is underdetermined, the minimal norm solution is returned. If A is rank deficient or singular then linear_solve returns []. If b is a matrix that is compatible with A then the problem is solved for the matrix valued right hand side and a matrix is returned. Note that if you want to solve Ax=b1 and Ax=b2 that you need to form the matrix transpose([b1,b2]) for the right hand side and then transpose the returned value. The solution is computed using QR factorization. If pivot is set to true (the default) then pivoting is used in the QR factorization, which is slower but expected to be more accurate.

Arguments:

Synopsis: Fast solution to Ax=b where A is 3x3.

Topics: Matrices, Linear Algebra

See Also: linear_solve(), matrix_inverse(), rot_inverse(), back_substitute(), cholesky()

Usage:

• x = linear_solve3(A,b)

Description:

Fast solution to a 3x3 linear system using Cramer's rule (which appears to be the fastest method in OpenSCAD). The input A must be a 3x3 matrix. Returns undef if A is singular. The input b must be a 3-vector. Note that Cramer's rule is not a stable algorithm, so for the highest accuracy on ill-conditioned problems you may want to use the general solver, which is about ten times slower.

Arguments:

Synopsis: General matrix inverse.

Topics: Matrices, Linear Algebra

See Also: linear_solve(), linear_solve3(), rot_inverse(), back_substitute(), cholesky()

Usage:

• mat = matrix_inverse(A)

Description:

Compute the matrix inverse of the square matrix A. If A is singular, returns undef. Note that if you just want to solve a linear system of equations you should NOT use this function. Instead use linear_solve(), or use qr_factor(). The computation will be faster and more accurate.

Synopsis: Invert 2d or 3d rotation transformations.

Topics: Matrices, Linear Algebra, Affine

See Also: linear_solve(), linear_solve3(), matrix_inverse(), back_substitute(), cholesky()

Usage:

• B = rot_inverse(A)

Description:

Inverts a 2d (3x3) or 3d (4x4) rotation matrix. The matrix can be a rotation around any center, so it may include a translation. This is faster and likely to be more accurate than using matrix_inverse().

Synopsis: Return basis for the null space of A.

Topics: Matrices, Linear Algebra

See Also: linear_solve(), linear_solve3(), matrix_inverse(), rot_inverse(), back_substitute(), cholesky()

Usage:

• x = null_space(A)

Description:

Returns an orthonormal basis for the null space of A, namely the vectors {x} such that Ax=0. If the null space is just the origin then returns an empty list.

Synopsis: Compute QR factorization of a matrix.

Topics: Matrices, Linear Algebra

See Also: linear_solve(), linear_solve3(), matrix_inverse(), rot_inverse(), back_substitute(), cholesky()

Usage:

• qr = qr_factor(A,[pivot]);

Description:

Calculates the QR factorization of the input matrix A and returns it as the list [Q,R,P]. This factorization can be used to solve linear systems of equations. The factorization is A = Q*R*transpose(P). If pivot is false (the default) then P is the identity matrix and A = Q*R. If pivot is true then column pivoting results in an R matrix where the diagonal is non-decreasing. The use of pivoting is supposed to increase accuracy for poorly conditioned problems, and is necessary for rank estimation or computation of the null space, but it may be slower.

Synopsis: Solve an upper triangular system, Rx=b.

Topics: Matrices, Linear Algebra

See Also: linear_solve(), linear_solve3(), matrix_inverse(), rot_inverse(), cholesky()

Usage:

• x = back_substitute(R, b, [transpose]);

Description:

Solves the problem Rx=b where R is an upper triangular square matrix. The lower triangular entries of R are ignored. If transpose==true then instead solve transpose(R)*x=b. You can supply a compatible matrix b and it will produce the solution for every column of b. Note that if you want to solve Rx=b1 and Rx=b2 you must set b to transpose([b1,b2]) and then take the transpose of the result. If the matrix is singular (e.g. has a zero on the diagonal) then it returns [].

Synopsis: Compute the Cholesky factorization of a matrix.

Topics: Matrices, Linear Algebra

See Also: linear_solve(), linear_solve3(), matrix_inverse(), rot_inverse(), back_substitute()

Usage:

• L = cholesky(A);

Description:

Compute the cholesky factor, L, of the symmetric positive definite matrix A. The matrix L is lower triangular and L * transpose(L) = A. If the A is not symmetric then an error is displayed. If the matrix is symmetric but not positive definite then undef is returned.

Synopsis: Compute determinant of 2x2 matrix.

Topics: Matrices, Linear Algebra

See Also: det3(), det4(), determinant(), norm_fro(), matrix_trace()

Usage:

• d = det2(M);

Description:

Rturns the determinant for the given 2x2 matrix.

Arguments:

Example 1:

include <BOSL2/std.scad>
M = [ [6,-2], [1,8] ];
det = det2(M);  // Returns: 50

Synopsis: Compute determinant of 3x3 matrix.

Topics: Matrices, Linear Algebra

See Also: det2(), det4(), determinant(), norm_fro(), matrix_trace()

Usage:

• d = det3(M);

Description:

Returns the determinant for the given 3x3 matrix.

Arguments:

Example 1:

include <BOSL2/std.scad>
M = [ [6,4,-2], [1,-2,8], [1,5,7] ];
det = det3(M);  // Returns: -334

Synopsis: Compute determinant of 4x4 matrix.

Topics: Matrices, Linear Algebra

See Also: det2(), det3(), determinant(), norm_fro(), matrix_trace()

Usage:

• d = det4(M);

Description:

Returns the determinant for the given 4x4 matrix.

Arguments:

Example 1:

include <BOSL2/std.scad>
M = [ [6,4,-2,1], [1,-2,8,-3], [1,5,7,4], [2,3,4,7] ];
det = det4(M);  // Returns: -1773

Synopsis: compute determinant of an arbitrary square matrix.

Topics: Matrices, Linear Algebra

See Also: det2(), det3(), det4(), norm_fro(), matrix_trace()

Usage:

• d = determinant(M);

Description:

Returns the determinant for the given square matrix.

Arguments:

Example 1:

include <BOSL2/std.scad>
M = [ [6,4,-2,9], [1,-2,8,3], [1,5,7,6], [4,2,5,1] ];
det = determinant(M);  // Returns: 2267

Synopsis: Compute Frobenius norm of a matrix

Topics: Matrices, Linear Algebra

See Also: det2(), det3(), det4(), determinant(), matrix_trace()

Usage:

• norm_fro(A)

Description:

Computes frobenius norm of input matrix. The frobenius norm is the square root of the sum of the squares of all of the entries of the matrix. On vectors it is the same as the usual 2-norm. This is an easily computed norm that is convenient for comparing two matrices.

Synopsis: Compute the trace of a square matrix.

Topics: Matrices, Linear Algebra

See Also: det2(), det3(), det4(), determinant(), norm_fro()

Usage:

• matrix_trace(M)

Description:

Computes the trace of a square matrix, the sum of the entries on the diagonal.