mirror of
https://github.com/hanzoai/jube-fork.git
synced 2026-10-30 18:33:33 +00:00
richard-churchman
cc0a2bee21
Included open-source Accord.net library code natively to the solution given the deprecated state of that project. Accord.net will be supported inside Jube for own needs (such as upgrading to .Net 7+). Accord.net libraries modified to compile to .Net 6. Modified Accord.net libraries to better support Newtonsoft Json.net serialisation. Included annotations in the Accord.net libraries to house important dependency statistics. Serialised json can be used independently of Jube to recall a Neural Network. Removal of Binary Serializer totally as a precursor to the projects support of .Net 7+. Extensive bug fixing in crucial functions of Neural Network training and recall. While this part of the system remains in Alpha, it is substantially more fit for purpose.
2206 lines
81 KiB
C#
2206 lines
81 KiB
C#
// Accord Math Library
|
|
// The Accord.NET Framework
|
|
// http://accord-framework.net
|
|
//
|
|
// Copyright © César Souza, 2009-2017
|
|
// cesarsouza at gmail.com
|
|
//
|
|
// This library is free software; you can redistribute it and/or
|
|
// modify it under the terms of the GNU Lesser General Public
|
|
// License as published by the Free Software Foundation; either
|
|
// version 2.1 of the License, or (at your option) any later version.
|
|
//
|
|
// This library is distributed in the hope that it will be useful,
|
|
// but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
// Lesser General Public License for more details.
|
|
//
|
|
// You should have received a copy of the GNU Lesser General Public
|
|
// License along with this library; if not, write to the Free Software
|
|
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
|
|
//
|
|
|
|
// ======================================================================
|
|
// This code has been generated by a tool; do not edit manually. Instead,
|
|
// edit the T4 template Matrix.Elementwise.tt so this file can be regenerated.
|
|
// ======================================================================
|
|
|
|
namespace Accord.Math
|
|
{
|
|
using System;
|
|
using System.CodeDom.Compiler;
|
|
using Accord.Math;
|
|
using System.Runtime.CompilerServices;
|
|
using Accord.Math.Decompositions;
|
|
|
|
//[GeneratedCode("Accord.NET T4 Templates", "3.1")]
|
|
public static partial class Matrix
|
|
{
|
|
|
|
/// <summary>
|
|
/// Returns the solution matrix if the matrix is square or the least squares solution otherwise.
|
|
/// </summary>
|
|
///
|
|
/// <param name="matrix">The matrix for the linear problem.</param>
|
|
/// <param name="rightSide">The right side <c>b</c>.</param>
|
|
/// <param name="leastSquares">True to produce a solution even if the
|
|
/// <paramref name="matrix"/> is singular; false otherwise. Default is false.</param>
|
|
///
|
|
/// <remarks>
|
|
/// Please note that this does not check if the matrix is non-singular
|
|
/// before attempting to solve. If a least squares solution is desired
|
|
/// in case the matrix is singular, pass true to the <paramref name="leastSquares"/>
|
|
/// parameter when calling this function.
|
|
/// </remarks>
|
|
///
|
|
/// <example>
|
|
/// <code>
|
|
/// // Create a matrix. Please note that this matrix
|
|
/// // is singular (i.e. not invertible), so only a
|
|
/// // least squares solution would be feasible here.
|
|
///
|
|
/// Double[,] matrix =
|
|
/// {
|
|
/// { 1, 2, 3 },
|
|
/// { 4, 5, 6 },
|
|
/// { 7, 8, 9 },
|
|
/// };
|
|
///
|
|
/// // Define a right side matrix b:
|
|
/// Double[,] rightSide = { {1}, {2}, {3} };
|
|
///
|
|
/// // Solve the linear system Ax = b by finding x:
|
|
/// Double[,] x = Matrix.Solve(matrix, rightSide, leastSquares: true);
|
|
///
|
|
/// // The answer should be { {-1/18}, {2/18}, {5/18} }.
|
|
/// </code>
|
|
/// </example>
|
|
///
|
|
public static Double[,] Solve(this Double[,] matrix, Double[,] rightSide, bool leastSquares = false)
|
|
{
|
|
int rows = matrix.GetLength(0);
|
|
int cols = matrix.GetLength(1);
|
|
|
|
if (rows != rightSide.GetLength(0))
|
|
throw new DimensionMismatchException("rightSide",
|
|
"The number of rows in the right hand side matrix must be "
|
|
+ "equal to the number of rows in the problem matrix.");
|
|
|
|
if (leastSquares)
|
|
{
|
|
return new SingularValueDecomposition(matrix,
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true).Solve(rightSide);
|
|
}
|
|
|
|
|
|
if (rows == cols)
|
|
{
|
|
// Solve by LU Decomposition if matrix is square.
|
|
return new LuDecomposition(matrix).Solve(rightSide);
|
|
}
|
|
else
|
|
{
|
|
if (cols < rows)
|
|
{
|
|
// Solve by QR Decomposition if not.
|
|
return new QrDecomposition(matrix).Solve(rightSide);
|
|
}
|
|
else
|
|
{
|
|
return new SingularValueDecomposition(matrix,
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true).Solve(rightSide);
|
|
}
|
|
}
|
|
}
|
|
|
|
/// <summary>
|
|
/// Returns the solution matrix if the matrix is square or the least squares solution otherwise.
|
|
/// </summary>
|
|
///
|
|
/// <param name="matrix">The matrix for the linear problem.</param>
|
|
/// <param name="rightSide">The right side <c>b</c>.</param>
|
|
/// <param name="leastSquares">True to produce a solution even if the
|
|
/// <paramref name="matrix"/> is singular; false otherwise. Default is false.</param>
|
|
///
|
|
/// <remarks>
|
|
/// Please note that this does not check if the matrix is non-singular
|
|
/// before attempting to solve. If a least squares solution is desired
|
|
/// in case the matrix is singular, pass true to the <paramref name="leastSquares"/>
|
|
/// parameter when calling this function.
|
|
/// </remarks>
|
|
///
|
|
/// <example>
|
|
/// <code>
|
|
/// // Create a matrix. Please note that this matrix
|
|
/// // is singular (i.e. not invertible), so only a
|
|
/// // least squares solution would be feasible here.
|
|
///
|
|
/// Double[,] matrix =
|
|
/// {
|
|
/// { 1, 2, 3 },
|
|
/// { 4, 5, 6 },
|
|
/// { 7, 8, 9 },
|
|
/// };
|
|
///
|
|
/// // Define a right side vector b:
|
|
/// Double[] rightSide = { 1, 2, 3 };
|
|
///
|
|
/// // Solve the linear system Ax = b by finding x:
|
|
/// Double[] x = Matrix.Solve(matrix, rightSide, leastSquares: true);
|
|
///
|
|
/// // The answer should be { -1/18, 2/18, 5/18 }.
|
|
/// </code>
|
|
/// </example>
|
|
///
|
|
public static Double[] Solve(this Double[,] matrix, Double[] rightSide, bool leastSquares = false)
|
|
{
|
|
if (matrix == null)
|
|
throw new ArgumentNullException("matrix");
|
|
|
|
if (rightSide == null)
|
|
throw new ArgumentNullException("rightSide");
|
|
|
|
int rows = matrix.GetLength(0);
|
|
int cols = matrix.GetLength(1);
|
|
|
|
if (rows != rightSide.Length)
|
|
throw new DimensionMismatchException("rightSide",
|
|
"The right hand side vector must have the same length"
|
|
+ "as there are rows of the problem matrix.");
|
|
|
|
if (leastSquares)
|
|
{
|
|
return new SingularValueDecomposition(matrix,
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true).Solve(rightSide);
|
|
}
|
|
|
|
|
|
if (rows == cols)
|
|
{
|
|
// Solve by LU Decomposition if matrix is square.
|
|
return new LuDecomposition(matrix).Solve(rightSide);
|
|
}
|
|
else
|
|
{
|
|
if (cols < rows)
|
|
{
|
|
// Solve by QR Decomposition if not.
|
|
return new QrDecomposition(matrix).Solve(rightSide);
|
|
}
|
|
else
|
|
{
|
|
return new SingularValueDecomposition(matrix,
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true).Solve(rightSide);
|
|
}
|
|
}
|
|
}
|
|
|
|
/// <summary>
|
|
/// Computes the inverse of a matrix.
|
|
/// </summary>
|
|
///
|
|
/// <example>
|
|
/// <code source="Unit Tests\Accord.Tests.Math\Matrix\MatrixTest.cs" region="doc_inverse" />
|
|
/// </example>
|
|
///
|
|
public static Double[,] Inverse(this Double[,] matrix)
|
|
{
|
|
return Inverse(matrix, false);
|
|
}
|
|
|
|
/// <summary>
|
|
/// Computes the inverse of a matrix.
|
|
/// </summary>
|
|
///
|
|
/// <example>
|
|
/// <code source="Unit Tests\Accord.Tests.Math\Matrix\MatrixTest.cs" region="doc_inverse" />
|
|
/// </example>
|
|
///
|
|
public static Double[,] Inverse(this Double[,] matrix, bool inPlace)
|
|
{
|
|
int rows = matrix.GetLength(0);
|
|
int cols = matrix.GetLength(1);
|
|
|
|
if (rows != cols)
|
|
throw new ArgumentException("Matrix must be square", "matrix");
|
|
|
|
if (rows == 3)
|
|
{
|
|
// Special case for 3x3 matrices
|
|
Double a = matrix[0, 0], b = matrix[0, 1], c = matrix[0, 2];
|
|
Double d = matrix[1, 0], e = matrix[1, 1], f = matrix[1, 2];
|
|
Double g = matrix[2, 0], h = matrix[2, 1], i = matrix[2, 2];
|
|
|
|
Double den = a * (e * i - f * h) -
|
|
b * (d * i - f * g) +
|
|
c * (d * h - e * g);
|
|
|
|
if (den == 0)
|
|
throw new SingularMatrixException();
|
|
|
|
Double m = 1 / den;
|
|
|
|
var inv = (inPlace) ? matrix : new Double[3, 3];
|
|
inv[0, 0] = m * (e * i - f * h);
|
|
inv[0, 1] = m * (c * h - b * i);
|
|
inv[0, 2] = m * (b * f - c * e);
|
|
inv[1, 0] = m * (f * g - d * i);
|
|
inv[1, 1] = m * (a * i - c * g);
|
|
inv[1, 2] = m * (c * d - a * f);
|
|
inv[2, 0] = m * (d * h - e * g);
|
|
inv[2, 1] = m * (b * g - a * h);
|
|
inv[2, 2] = m * (a * e - b * d);
|
|
|
|
return inv;
|
|
}
|
|
|
|
if (rows == 2)
|
|
{
|
|
// Special case for 2x2 matrices
|
|
Double a = matrix[0, 0], b = matrix[0, 1];
|
|
Double c = matrix[1, 0], d = matrix[1, 1];
|
|
|
|
Double den = a * d - b * c;
|
|
|
|
if (den == 0)
|
|
throw new SingularMatrixException();
|
|
|
|
Double m = 1 / den;
|
|
|
|
var inv = (inPlace) ? matrix : new Double[2, 2];
|
|
inv[0, 0] = +m * d;
|
|
inv[0, 1] = -m * b;
|
|
inv[1, 0] = -m * c;
|
|
inv[1, 1] = +m * a;
|
|
|
|
return inv;
|
|
}
|
|
|
|
return new LuDecomposition(matrix, false, inPlace).Inverse();
|
|
}
|
|
|
|
/// <summary>
|
|
/// Computes the pseudo-inverse of a matrix.
|
|
/// </summary>
|
|
///
|
|
/// <example>
|
|
/// <code source="Unit Tests\Accord.Tests.Math\Matrix\MatrixTest.cs" region="doc_pseudoinverse" />
|
|
/// </example>
|
|
///
|
|
public static Double[,] PseudoInverse(this Double[,] matrix)
|
|
{
|
|
return new SingularValueDecomposition(matrix,
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true).Inverse();
|
|
}
|
|
|
|
|
|
/// <summary>
|
|
/// Returns the solution matrix if the matrix is square or the least squares solution otherwise.
|
|
/// </summary>
|
|
///
|
|
/// <param name="matrix">The matrix for the linear problem.</param>
|
|
/// <param name="rightSide">The right side <c>b</c>.</param>
|
|
/// <param name="leastSquares">True to produce a solution even if the
|
|
/// <paramref name="matrix"/> is singular; false otherwise. Default is false.</param>
|
|
///
|
|
/// <remarks>
|
|
/// Please note that this does not check if the matrix is non-singular
|
|
/// before attempting to solve. If a least squares solution is desired
|
|
/// in case the matrix is singular, pass true to the <paramref name="leastSquares"/>
|
|
/// parameter when calling this function.
|
|
/// </remarks>
|
|
///
|
|
/// <example>
|
|
/// <code>
|
|
/// // Create a matrix. Please note that this matrix
|
|
/// // is singular (i.e. not invertible), so only a
|
|
/// // least squares solution would be feasible here.
|
|
///
|
|
/// Double[,] matrix =
|
|
/// {
|
|
/// { 1, 2, 3 },
|
|
/// { 4, 5, 6 },
|
|
/// { 7, 8, 9 },
|
|
/// };
|
|
///
|
|
/// // Define a right side matrix b:
|
|
/// Double[,] rightSide = { {1}, {2}, {3} };
|
|
///
|
|
/// // Solve the linear system Ax = b by finding x:
|
|
/// Double[,] x = Matrix.Solve(matrix, rightSide, leastSquares: true);
|
|
///
|
|
/// // The answer should be { {-1/18}, {2/18}, {5/18} }.
|
|
/// </code>
|
|
/// </example>
|
|
///
|
|
public static Double[][] Solve(this Double[][] matrix, Double[][] rightSide, bool leastSquares = false)
|
|
{
|
|
if (matrix.Length != rightSide.Length)
|
|
{
|
|
throw new DimensionMismatchException("rightSide",
|
|
"The number of rows in the right hand side matrix must be "
|
|
+ "equal to the number of rows in the problem matrix.");
|
|
}
|
|
|
|
return matrix.Decompose(leastSquares).Solve(rightSide);
|
|
}
|
|
|
|
/// <summary>
|
|
/// Returns the solution matrix if the matrix is square or the least squares solution otherwise.
|
|
/// </summary>
|
|
///
|
|
/// <param name="matrix">The matrix for the linear problem.</param>
|
|
/// <param name="rightSide">The right side <c>b</c>.</param>
|
|
/// <param name="leastSquares">True to produce a solution even if the
|
|
/// <paramref name="matrix"/> is singular; false otherwise. Default is false.</param>
|
|
///
|
|
/// <remarks>
|
|
/// Please note that this does not check if the matrix is non-singular
|
|
/// before attempting to solve. If a least squares solution is desired
|
|
/// in case the matrix is singular, pass true to the <paramref name="leastSquares"/>
|
|
/// parameter when calling this function.
|
|
/// </remarks>
|
|
///
|
|
/// <example>
|
|
/// <code>
|
|
/// // Create a matrix. Please note that this matrix
|
|
/// // is singular (i.e. not invertible), so only a
|
|
/// // least squares solution would be feasible here.
|
|
///
|
|
/// Double[,] matrix =
|
|
/// {
|
|
/// { 1, 2, 3 },
|
|
/// { 4, 5, 6 },
|
|
/// { 7, 8, 9 },
|
|
/// };
|
|
///
|
|
/// // Define a right side vector b:
|
|
/// Double[] rightSide = { 1, 2, 3 };
|
|
///
|
|
/// // Solve the linear system Ax = b by finding x:
|
|
/// Double[] x = Matrix.Solve(matrix, rightSide, leastSquares: true);
|
|
///
|
|
/// // The answer should be { -1/18, 2/18, 5/18 }.
|
|
/// </code>
|
|
/// </example>
|
|
///
|
|
public static Double[] Solve(this Double[][] matrix, Double[] rightSide, bool leastSquares = false)
|
|
{
|
|
if (matrix.Length != rightSide.Length)
|
|
{
|
|
throw new DimensionMismatchException("rightSide",
|
|
"The right hand side vector must have the same length"
|
|
+ "as there are rows of the problem matrix.");
|
|
}
|
|
|
|
return matrix.Decompose(leastSquares).Solve(rightSide);
|
|
}
|
|
|
|
/// <summary>
|
|
/// Returns the solution matrix for a linear system involving a diagonal matrix ion the right-hand side.
|
|
/// </summary>
|
|
///
|
|
/// <param name="matrix">The matrix for the linear problem.</param>
|
|
/// <param name="diagonalRightSide">The right side <c>b</c>.</param>
|
|
/// <param name="leastSquares">True to produce a solution even if the
|
|
/// <paramref name="matrix"/> is singular; false otherwise. Default is false.</param>
|
|
///
|
|
/// <remarks>
|
|
/// Please note that this does not check if the matrix is non-singular
|
|
/// before attempting to solve. If a least squares solution is desired
|
|
/// in case the matrix is singular, pass true to the <paramref name="leastSquares"/>
|
|
/// parameter when calling this function.
|
|
/// </remarks>
|
|
///
|
|
public static Double[][] SolveForDiagonal(this Double[][] matrix, Double[] diagonalRightSide, bool leastSquares = false)
|
|
{
|
|
if (matrix.Length != diagonalRightSide.Length)
|
|
{
|
|
throw new DimensionMismatchException("diagonalRightSide",
|
|
"The right hand side matrix must have the same length"
|
|
+ "as there are rows of the problem matrix.");
|
|
}
|
|
|
|
return matrix.Decompose(leastSquares).SolveForDiagonal(diagonalRightSide);
|
|
}
|
|
|
|
|
|
/// <summary>
|
|
/// Creates a matrix decomposition that be used to compute the solution matrix if the
|
|
/// matrix is square or the least squares solution otherwise.
|
|
/// </summary>
|
|
///
|
|
public static ISolverMatrixDecomposition<Double> Decompose(this Double[,] matrix, bool leastSquares = false)
|
|
{
|
|
int rows = matrix.Rows();
|
|
int cols = matrix.Columns();
|
|
|
|
if (leastSquares)
|
|
{
|
|
return new SingularValueDecomposition(matrix,
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true);
|
|
}
|
|
|
|
if (rows == cols)
|
|
{
|
|
// Solve by LU Decomposition if matrix is square.
|
|
return new LuDecomposition(matrix);
|
|
}
|
|
else
|
|
{
|
|
if (cols < rows)
|
|
{
|
|
// Solve by QR Decomposition if not.
|
|
return new QrDecomposition(matrix);
|
|
}
|
|
else
|
|
{
|
|
return new SingularValueDecomposition(matrix,
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true);
|
|
}
|
|
}
|
|
}
|
|
|
|
/// <summary>
|
|
/// Creates a matrix decomposition that be used to compute the solution matrix if the
|
|
/// matrix is square or the least squares solution otherwise.
|
|
/// </summary>
|
|
///
|
|
public static ISolverArrayDecomposition<Double> Decompose(this Double[][] matrix, bool leastSquares = false)
|
|
{
|
|
int rows = matrix.Rows();
|
|
int cols = matrix.Columns();
|
|
|
|
if (leastSquares)
|
|
{
|
|
return new JaggedSingularValueDecomposition(matrix,
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true);
|
|
}
|
|
|
|
if (rows == cols)
|
|
{
|
|
// Solve by LU Decomposition if matrix is square.
|
|
return new JaggedLuDecomposition(matrix);
|
|
}
|
|
else
|
|
{
|
|
if (cols < rows)
|
|
{
|
|
// Solve by QR Decomposition if not.
|
|
return new JaggedQrDecomposition(matrix);
|
|
}
|
|
else
|
|
{
|
|
return new JaggedSingularValueDecomposition(matrix,
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true);
|
|
}
|
|
}
|
|
}
|
|
|
|
/// <summary>
|
|
/// Computes the inverse of a matrix.
|
|
/// </summary>
|
|
///
|
|
///
|
|
/// <example>
|
|
/// <code source="Unit Tests\Accord.Tests.Math\Matrix\MatrixTest.cs" region="doc_inverse" />
|
|
/// </example>
|
|
///
|
|
public static Double[][] Inverse(this Double[][] matrix)
|
|
{
|
|
return Inverse(matrix, false);
|
|
}
|
|
|
|
/// <summary>
|
|
/// Computes the inverse of a matrix.
|
|
/// </summary>
|
|
///
|
|
/// <example>
|
|
/// <code source="Unit Tests\Accord.Tests.Math\Matrix\MatrixTest.cs" region="doc_inverse" />
|
|
/// </example>
|
|
///
|
|
public static Double[][] Inverse(this Double[][] matrix, bool inPlace)
|
|
{
|
|
int rows = matrix.Length;
|
|
int cols = matrix[0].Length;
|
|
|
|
if (rows != cols)
|
|
throw new ArgumentException("Matrix must be square", "matrix");
|
|
|
|
if (rows == 3)
|
|
{
|
|
// Special case for 3x3 matrices
|
|
Double a = matrix[0][0], b = matrix[0][1], c = matrix[0][2];
|
|
Double d = matrix[1][0], e = matrix[1][1], f = matrix[1][2];
|
|
Double g = matrix[2][0], h = matrix[2][1], i = matrix[2][2];
|
|
|
|
Double den = a * (e * i - f * h) -
|
|
b * (d * i - f * g) +
|
|
c * (d * h - e * g);
|
|
|
|
if (den == 0)
|
|
throw new SingularMatrixException();
|
|
|
|
Double m = 1 / den;
|
|
|
|
var inv = matrix;
|
|
if (!inPlace)
|
|
{
|
|
inv = new Double[3][];
|
|
for (int j = 0; j < inv.Length; j++)
|
|
inv[j] = new Double[3];
|
|
}
|
|
|
|
inv[0][0] = m * (e * i - f * h);
|
|
inv[0][1] = m * (c * h - b * i);
|
|
inv[0][2] = m * (b * f - c * e);
|
|
inv[1][0] = m * (f * g - d * i);
|
|
inv[1][1] = m * (a * i - c * g);
|
|
inv[1][2] = m * (c * d - a * f);
|
|
inv[2][0] = m * (d * h - e * g);
|
|
inv[2][1] = m * (b * g - a * h);
|
|
inv[2][2] = m * (a * e - b * d);
|
|
|
|
return inv;
|
|
}
|
|
|
|
if (rows == 2)
|
|
{
|
|
// Special case for 2x2 matrices
|
|
Double a = matrix[0][0], b = matrix[0][1];
|
|
Double c = matrix[1][0], d = matrix[1][1];
|
|
|
|
Double den = a * d - b * c;
|
|
|
|
if (den == 0)
|
|
throw new SingularMatrixException();
|
|
|
|
Double m = 1 / den;
|
|
|
|
var inv = matrix;
|
|
if (!inPlace)
|
|
{
|
|
inv = new Double[2][];
|
|
for (int j = 0; j < inv.Length; j++)
|
|
inv[j] = new Double[2];
|
|
}
|
|
|
|
inv[0][0] = +m * d;
|
|
inv[0][1] = -m * b;
|
|
inv[1][0] = -m * c;
|
|
inv[1][1] = +m * a;
|
|
|
|
return inv;
|
|
}
|
|
|
|
return new JaggedLuDecomposition(matrix, false, inPlace).Inverse();
|
|
}
|
|
|
|
/// <summary>
|
|
/// Computes the pseudo-inverse of a matrix.
|
|
/// </summary>
|
|
///
|
|
/// <example>
|
|
/// <code source="Unit Tests\Accord.Tests.Math\Matrix\MatrixTest.cs" region="doc_pseudoinverse" />
|
|
/// </example>
|
|
///
|
|
public static Double[][] PseudoInverse(this Double[][] matrix)
|
|
{
|
|
return new JaggedSingularValueDecomposition(matrix,
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true).Inverse();
|
|
}
|
|
|
|
/// <summary>
|
|
/// Divides two matrices by multiplying A by the inverse of B.
|
|
/// </summary>
|
|
///
|
|
/// <param name="a">The first matrix.</param>
|
|
/// <param name="b">The second matrix (which will be inverted).</param>
|
|
/// <param name="leastSquares">True to produce a solution even if the
|
|
/// <paramref name="b"/> is singular; false otherwise. Default is false.</param>
|
|
///
|
|
/// <returns>The result from the division <c>AB^-1</c> of the given matrices.</returns>
|
|
///
|
|
public static Double[,] Divide(this Double[,] a, Double[,] b, bool leastSquares = false)
|
|
{
|
|
int rows = b.Rows();
|
|
int cols = b.Columns();
|
|
|
|
if (leastSquares)
|
|
{
|
|
return new SingularValueDecomposition(b.Transpose(),
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true).Solve(a.Transpose()).Transpose();
|
|
}
|
|
|
|
|
|
if (rows == cols && cols == a.Rows())
|
|
{
|
|
// Solve by LU Decomposition if matrix is square.
|
|
return new LuDecomposition(b, true).SolveTranspose(a);
|
|
}
|
|
else
|
|
{
|
|
if (rows <= cols)
|
|
{
|
|
// Solve by QR Decomposition if not.
|
|
return new QrDecomposition(b, true).SolveTranspose(a);
|
|
}
|
|
else
|
|
{
|
|
return new SingularValueDecomposition(b.Transpose(),
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true).Solve(a.Transpose()).Transpose();
|
|
}
|
|
}
|
|
}
|
|
|
|
/// <summary>
|
|
/// Divides two matrices by multiplying A by the inverse of B.
|
|
/// </summary>
|
|
///
|
|
/// <param name="a">The first matrix.</param>
|
|
/// <param name="b">The second matrix (which will be inverted).</param>
|
|
/// <param name="leastSquares">True to produce a solution even if the
|
|
/// <paramref name="b"/> is singular; false otherwise. Default is false.</param>
|
|
///
|
|
/// <returns>The result from the division <c>AB^-1</c> of the given matrices.</returns>
|
|
///
|
|
public static Double[][] Divide(this Double[][] a, Double[][] b, bool leastSquares = false)
|
|
{
|
|
int rows = b.Rows();
|
|
int cols = b.Columns();
|
|
|
|
if (leastSquares)
|
|
{
|
|
return new JaggedSingularValueDecomposition(b.Transpose(),
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true).Solve(a.Transpose()).Transpose();
|
|
}
|
|
|
|
|
|
if (rows == cols && cols == a.Rows())
|
|
{
|
|
// Solve by LU Decomposition if matrix is square.
|
|
return new JaggedLuDecomposition(b, true).SolveTranspose(a);
|
|
}
|
|
else
|
|
{
|
|
if (rows <= cols)
|
|
{
|
|
// Solve by QR Decomposition if not.
|
|
return new JaggedQrDecomposition(b, true).SolveTranspose(a);
|
|
}
|
|
else
|
|
{
|
|
return new JaggedSingularValueDecomposition(b.Transpose(),
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true).Solve(a.Transpose()).Transpose();
|
|
}
|
|
}
|
|
}
|
|
|
|
/// <summary>
|
|
/// Gets the null-space of a column vector.
|
|
/// </summary>
|
|
///
|
|
public static Double[][] Null(this Double[] vector)
|
|
{
|
|
return Null(Jagged.ColumnVector(vector));
|
|
}
|
|
|
|
/// <summary>
|
|
/// Gets the null-space of a matrix.
|
|
/// </summary>
|
|
///
|
|
public static Double[][] Null(this Double[][] matrix)
|
|
{
|
|
var qr = new JaggedQrDecomposition(matrix, economy: false);
|
|
var Q = qr.OrthogonalFactor;
|
|
var threshold = matrix.GetLength().Max() * Constants.DoubleEpsilon;
|
|
int[] idx = qr.Diagonal.Find(x => (Double)Math.Abs(x) < threshold);
|
|
return Q.GetColumns(idx);
|
|
}
|
|
|
|
/// <summary>
|
|
/// Gets the null-space of a matrix.
|
|
/// </summary>
|
|
///
|
|
public static Double[,] Null(this Double[,] matrix)
|
|
{
|
|
var qr = new QrDecomposition(matrix, economy: false);
|
|
var Q = qr.OrthogonalFactor;
|
|
var threshold = matrix.GetLength().Max() * Constants.DoubleEpsilon;
|
|
int[] idx = qr.Diagonal.Find(x => (Double)Math.Abs(x) < threshold);
|
|
return Q.GetColumns(idx);
|
|
}
|
|
|
|
/// <summary>
|
|
/// Returns the solution matrix if the matrix is square or the least squares solution otherwise.
|
|
/// </summary>
|
|
///
|
|
/// <param name="matrix">The matrix for the linear problem.</param>
|
|
/// <param name="rightSide">The right side <c>b</c>.</param>
|
|
/// <param name="leastSquares">True to produce a solution even if the
|
|
/// <paramref name="matrix"/> is singular; false otherwise. Default is false.</param>
|
|
///
|
|
/// <remarks>
|
|
/// Please note that this does not check if the matrix is non-singular
|
|
/// before attempting to solve. If a least squares solution is desired
|
|
/// in case the matrix is singular, pass true to the <paramref name="leastSquares"/>
|
|
/// parameter when calling this function.
|
|
/// </remarks>
|
|
///
|
|
/// <example>
|
|
/// <code>
|
|
/// // Create a matrix. Please note that this matrix
|
|
/// // is singular (i.e. not invertible), so only a
|
|
/// // least squares solution would be feasible here.
|
|
///
|
|
/// Single[,] matrix =
|
|
/// {
|
|
/// { 1, 2, 3 },
|
|
/// { 4, 5, 6 },
|
|
/// { 7, 8, 9 },
|
|
/// };
|
|
///
|
|
/// // Define a right side matrix b:
|
|
/// Single[,] rightSide = { {1}, {2}, {3} };
|
|
///
|
|
/// // Solve the linear system Ax = b by finding x:
|
|
/// Single[,] x = Matrix.Solve(matrix, rightSide, leastSquares: true);
|
|
///
|
|
/// // The answer should be { {-1/18}, {2/18}, {5/18} }.
|
|
/// </code>
|
|
/// </example>
|
|
///
|
|
public static Single[,] Solve(this Single[,] matrix, Single[,] rightSide, bool leastSquares = false)
|
|
{
|
|
int rows = matrix.GetLength(0);
|
|
int cols = matrix.GetLength(1);
|
|
|
|
if (rows != rightSide.GetLength(0))
|
|
throw new DimensionMismatchException("rightSide",
|
|
"The number of rows in the right hand side matrix must be "
|
|
+ "equal to the number of rows in the problem matrix.");
|
|
|
|
if (leastSquares)
|
|
{
|
|
return new SingularValueDecompositionF(matrix,
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true).Solve(rightSide);
|
|
}
|
|
|
|
|
|
if (rows == cols)
|
|
{
|
|
// Solve by LU Decomposition if matrix is square.
|
|
return new LuDecompositionF(matrix).Solve(rightSide);
|
|
}
|
|
else
|
|
{
|
|
if (cols < rows)
|
|
{
|
|
// Solve by QR Decomposition if not.
|
|
return new QrDecompositionF(matrix).Solve(rightSide);
|
|
}
|
|
else
|
|
{
|
|
return new SingularValueDecompositionF(matrix,
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true).Solve(rightSide);
|
|
}
|
|
}
|
|
}
|
|
|
|
/// <summary>
|
|
/// Returns the solution matrix if the matrix is square or the least squares solution otherwise.
|
|
/// </summary>
|
|
///
|
|
/// <param name="matrix">The matrix for the linear problem.</param>
|
|
/// <param name="rightSide">The right side <c>b</c>.</param>
|
|
/// <param name="leastSquares">True to produce a solution even if the
|
|
/// <paramref name="matrix"/> is singular; false otherwise. Default is false.</param>
|
|
///
|
|
/// <remarks>
|
|
/// Please note that this does not check if the matrix is non-singular
|
|
/// before attempting to solve. If a least squares solution is desired
|
|
/// in case the matrix is singular, pass true to the <paramref name="leastSquares"/>
|
|
/// parameter when calling this function.
|
|
/// </remarks>
|
|
///
|
|
/// <example>
|
|
/// <code>
|
|
/// // Create a matrix. Please note that this matrix
|
|
/// // is singular (i.e. not invertible), so only a
|
|
/// // least squares solution would be feasible here.
|
|
///
|
|
/// Single[,] matrix =
|
|
/// {
|
|
/// { 1, 2, 3 },
|
|
/// { 4, 5, 6 },
|
|
/// { 7, 8, 9 },
|
|
/// };
|
|
///
|
|
/// // Define a right side vector b:
|
|
/// Single[] rightSide = { 1, 2, 3 };
|
|
///
|
|
/// // Solve the linear system Ax = b by finding x:
|
|
/// Single[] x = Matrix.Solve(matrix, rightSide, leastSquares: true);
|
|
///
|
|
/// // The answer should be { -1/18, 2/18, 5/18 }.
|
|
/// </code>
|
|
/// </example>
|
|
///
|
|
public static Single[] Solve(this Single[,] matrix, Single[] rightSide, bool leastSquares = false)
|
|
{
|
|
if (matrix == null)
|
|
throw new ArgumentNullException("matrix");
|
|
|
|
if (rightSide == null)
|
|
throw new ArgumentNullException("rightSide");
|
|
|
|
int rows = matrix.GetLength(0);
|
|
int cols = matrix.GetLength(1);
|
|
|
|
if (rows != rightSide.Length)
|
|
throw new DimensionMismatchException("rightSide",
|
|
"The right hand side vector must have the same length"
|
|
+ "as there are rows of the problem matrix.");
|
|
|
|
if (leastSquares)
|
|
{
|
|
return new SingularValueDecompositionF(matrix,
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true).Solve(rightSide);
|
|
}
|
|
|
|
|
|
if (rows == cols)
|
|
{
|
|
// Solve by LU Decomposition if matrix is square.
|
|
return new LuDecompositionF(matrix).Solve(rightSide);
|
|
}
|
|
else
|
|
{
|
|
if (cols < rows)
|
|
{
|
|
// Solve by QR Decomposition if not.
|
|
return new QrDecompositionF(matrix).Solve(rightSide);
|
|
}
|
|
else
|
|
{
|
|
return new SingularValueDecompositionF(matrix,
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true).Solve(rightSide);
|
|
}
|
|
}
|
|
}
|
|
|
|
/// <summary>
|
|
/// Computes the inverse of a matrix.
|
|
/// </summary>
|
|
///
|
|
/// <example>
|
|
/// <code source="Unit Tests\Accord.Tests.Math\Matrix\MatrixTest.cs" region="doc_inverse" />
|
|
/// </example>
|
|
///
|
|
public static Single[,] Inverse(this Single[,] matrix)
|
|
{
|
|
return Inverse(matrix, false);
|
|
}
|
|
|
|
/// <summary>
|
|
/// Computes the inverse of a matrix.
|
|
/// </summary>
|
|
///
|
|
/// <example>
|
|
/// <code source="Unit Tests\Accord.Tests.Math\Matrix\MatrixTest.cs" region="doc_inverse" />
|
|
/// </example>
|
|
///
|
|
public static Single[,] Inverse(this Single[,] matrix, bool inPlace)
|
|
{
|
|
int rows = matrix.GetLength(0);
|
|
int cols = matrix.GetLength(1);
|
|
|
|
if (rows != cols)
|
|
throw new ArgumentException("Matrix must be square", "matrix");
|
|
|
|
if (rows == 3)
|
|
{
|
|
// Special case for 3x3 matrices
|
|
Single a = matrix[0, 0], b = matrix[0, 1], c = matrix[0, 2];
|
|
Single d = matrix[1, 0], e = matrix[1, 1], f = matrix[1, 2];
|
|
Single g = matrix[2, 0], h = matrix[2, 1], i = matrix[2, 2];
|
|
|
|
Single den = a * (e * i - f * h) -
|
|
b * (d * i - f * g) +
|
|
c * (d * h - e * g);
|
|
|
|
if (den == 0)
|
|
throw new SingularMatrixException();
|
|
|
|
Single m = 1 / den;
|
|
|
|
var inv = (inPlace) ? matrix : new Single[3, 3];
|
|
inv[0, 0] = m * (e * i - f * h);
|
|
inv[0, 1] = m * (c * h - b * i);
|
|
inv[0, 2] = m * (b * f - c * e);
|
|
inv[1, 0] = m * (f * g - d * i);
|
|
inv[1, 1] = m * (a * i - c * g);
|
|
inv[1, 2] = m * (c * d - a * f);
|
|
inv[2, 0] = m * (d * h - e * g);
|
|
inv[2, 1] = m * (b * g - a * h);
|
|
inv[2, 2] = m * (a * e - b * d);
|
|
|
|
return inv;
|
|
}
|
|
|
|
if (rows == 2)
|
|
{
|
|
// Special case for 2x2 matrices
|
|
Single a = matrix[0, 0], b = matrix[0, 1];
|
|
Single c = matrix[1, 0], d = matrix[1, 1];
|
|
|
|
Single den = a * d - b * c;
|
|
|
|
if (den == 0)
|
|
throw new SingularMatrixException();
|
|
|
|
Single m = 1 / den;
|
|
|
|
var inv = (inPlace) ? matrix : new Single[2, 2];
|
|
inv[0, 0] = +m * d;
|
|
inv[0, 1] = -m * b;
|
|
inv[1, 0] = -m * c;
|
|
inv[1, 1] = +m * a;
|
|
|
|
return inv;
|
|
}
|
|
|
|
return new LuDecompositionF(matrix, false, inPlace).Inverse();
|
|
}
|
|
|
|
/// <summary>
|
|
/// Computes the pseudo-inverse of a matrix.
|
|
/// </summary>
|
|
///
|
|
/// <example>
|
|
/// <code source="Unit Tests\Accord.Tests.Math\Matrix\MatrixTest.cs" region="doc_pseudoinverse" />
|
|
/// </example>
|
|
///
|
|
public static Single[,] PseudoInverse(this Single[,] matrix)
|
|
{
|
|
return new SingularValueDecompositionF(matrix,
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true).Inverse();
|
|
}
|
|
|
|
|
|
/// <summary>
|
|
/// Returns the solution matrix if the matrix is square or the least squares solution otherwise.
|
|
/// </summary>
|
|
///
|
|
/// <param name="matrix">The matrix for the linear problem.</param>
|
|
/// <param name="rightSide">The right side <c>b</c>.</param>
|
|
/// <param name="leastSquares">True to produce a solution even if the
|
|
/// <paramref name="matrix"/> is singular; false otherwise. Default is false.</param>
|
|
///
|
|
/// <remarks>
|
|
/// Please note that this does not check if the matrix is non-singular
|
|
/// before attempting to solve. If a least squares solution is desired
|
|
/// in case the matrix is singular, pass true to the <paramref name="leastSquares"/>
|
|
/// parameter when calling this function.
|
|
/// </remarks>
|
|
///
|
|
/// <example>
|
|
/// <code>
|
|
/// // Create a matrix. Please note that this matrix
|
|
/// // is singular (i.e. not invertible), so only a
|
|
/// // least squares solution would be feasible here.
|
|
///
|
|
/// Single[,] matrix =
|
|
/// {
|
|
/// { 1, 2, 3 },
|
|
/// { 4, 5, 6 },
|
|
/// { 7, 8, 9 },
|
|
/// };
|
|
///
|
|
/// // Define a right side matrix b:
|
|
/// Single[,] rightSide = { {1}, {2}, {3} };
|
|
///
|
|
/// // Solve the linear system Ax = b by finding x:
|
|
/// Single[,] x = Matrix.Solve(matrix, rightSide, leastSquares: true);
|
|
///
|
|
/// // The answer should be { {-1/18}, {2/18}, {5/18} }.
|
|
/// </code>
|
|
/// </example>
|
|
///
|
|
public static Single[][] Solve(this Single[][] matrix, Single[][] rightSide, bool leastSquares = false)
|
|
{
|
|
if (matrix.Length != rightSide.Length)
|
|
{
|
|
throw new DimensionMismatchException("rightSide",
|
|
"The number of rows in the right hand side matrix must be "
|
|
+ "equal to the number of rows in the problem matrix.");
|
|
}
|
|
|
|
return matrix.Decompose(leastSquares).Solve(rightSide);
|
|
}
|
|
|
|
/// <summary>
|
|
/// Returns the solution matrix if the matrix is square or the least squares solution otherwise.
|
|
/// </summary>
|
|
///
|
|
/// <param name="matrix">The matrix for the linear problem.</param>
|
|
/// <param name="rightSide">The right side <c>b</c>.</param>
|
|
/// <param name="leastSquares">True to produce a solution even if the
|
|
/// <paramref name="matrix"/> is singular; false otherwise. Default is false.</param>
|
|
///
|
|
/// <remarks>
|
|
/// Please note that this does not check if the matrix is non-singular
|
|
/// before attempting to solve. If a least squares solution is desired
|
|
/// in case the matrix is singular, pass true to the <paramref name="leastSquares"/>
|
|
/// parameter when calling this function.
|
|
/// </remarks>
|
|
///
|
|
/// <example>
|
|
/// <code>
|
|
/// // Create a matrix. Please note that this matrix
|
|
/// // is singular (i.e. not invertible), so only a
|
|
/// // least squares solution would be feasible here.
|
|
///
|
|
/// Single[,] matrix =
|
|
/// {
|
|
/// { 1, 2, 3 },
|
|
/// { 4, 5, 6 },
|
|
/// { 7, 8, 9 },
|
|
/// };
|
|
///
|
|
/// // Define a right side vector b:
|
|
/// Single[] rightSide = { 1, 2, 3 };
|
|
///
|
|
/// // Solve the linear system Ax = b by finding x:
|
|
/// Single[] x = Matrix.Solve(matrix, rightSide, leastSquares: true);
|
|
///
|
|
/// // The answer should be { -1/18, 2/18, 5/18 }.
|
|
/// </code>
|
|
/// </example>
|
|
///
|
|
public static Single[] Solve(this Single[][] matrix, Single[] rightSide, bool leastSquares = false)
|
|
{
|
|
if (matrix.Length != rightSide.Length)
|
|
{
|
|
throw new DimensionMismatchException("rightSide",
|
|
"The right hand side vector must have the same length"
|
|
+ "as there are rows of the problem matrix.");
|
|
}
|
|
|
|
return matrix.Decompose(leastSquares).Solve(rightSide);
|
|
}
|
|
|
|
/// <summary>
|
|
/// Returns the solution matrix for a linear system involving a diagonal matrix ion the right-hand side.
|
|
/// </summary>
|
|
///
|
|
/// <param name="matrix">The matrix for the linear problem.</param>
|
|
/// <param name="diagonalRightSide">The right side <c>b</c>.</param>
|
|
/// <param name="leastSquares">True to produce a solution even if the
|
|
/// <paramref name="matrix"/> is singular; false otherwise. Default is false.</param>
|
|
///
|
|
/// <remarks>
|
|
/// Please note that this does not check if the matrix is non-singular
|
|
/// before attempting to solve. If a least squares solution is desired
|
|
/// in case the matrix is singular, pass true to the <paramref name="leastSquares"/>
|
|
/// parameter when calling this function.
|
|
/// </remarks>
|
|
///
|
|
public static Single[][] SolveForDiagonal(this Single[][] matrix, Single[] diagonalRightSide, bool leastSquares = false)
|
|
{
|
|
if (matrix.Length != diagonalRightSide.Length)
|
|
{
|
|
throw new DimensionMismatchException("diagonalRightSide",
|
|
"The right hand side matrix must have the same length"
|
|
+ "as there are rows of the problem matrix.");
|
|
}
|
|
|
|
return matrix.Decompose(leastSquares).SolveForDiagonal(diagonalRightSide);
|
|
}
|
|
|
|
|
|
/// <summary>
|
|
/// Creates a matrix decomposition that be used to compute the solution matrix if the
|
|
/// matrix is square or the least squares solution otherwise.
|
|
/// </summary>
|
|
///
|
|
public static ISolverMatrixDecomposition<Single> Decompose(this Single[,] matrix, bool leastSquares = false)
|
|
{
|
|
int rows = matrix.Rows();
|
|
int cols = matrix.Columns();
|
|
|
|
if (leastSquares)
|
|
{
|
|
return new SingularValueDecompositionF(matrix,
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true);
|
|
}
|
|
|
|
if (rows == cols)
|
|
{
|
|
// Solve by LU Decomposition if matrix is square.
|
|
return new LuDecompositionF(matrix);
|
|
}
|
|
else
|
|
{
|
|
if (cols < rows)
|
|
{
|
|
// Solve by QR Decomposition if not.
|
|
return new QrDecompositionF(matrix);
|
|
}
|
|
else
|
|
{
|
|
return new SingularValueDecompositionF(matrix,
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true);
|
|
}
|
|
}
|
|
}
|
|
|
|
/// <summary>
|
|
/// Creates a matrix decomposition that be used to compute the solution matrix if the
|
|
/// matrix is square or the least squares solution otherwise.
|
|
/// </summary>
|
|
///
|
|
public static ISolverArrayDecomposition<Single> Decompose(this Single[][] matrix, bool leastSquares = false)
|
|
{
|
|
int rows = matrix.Rows();
|
|
int cols = matrix.Columns();
|
|
|
|
if (leastSquares)
|
|
{
|
|
return new JaggedSingularValueDecompositionF(matrix,
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true);
|
|
}
|
|
|
|
if (rows == cols)
|
|
{
|
|
// Solve by LU Decomposition if matrix is square.
|
|
return new JaggedLuDecompositionF(matrix);
|
|
}
|
|
else
|
|
{
|
|
if (cols < rows)
|
|
{
|
|
// Solve by QR Decomposition if not.
|
|
return new JaggedQrDecompositionF(matrix);
|
|
}
|
|
else
|
|
{
|
|
return new JaggedSingularValueDecompositionF(matrix,
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true);
|
|
}
|
|
}
|
|
}
|
|
|
|
/// <summary>
|
|
/// Computes the inverse of a matrix.
|
|
/// </summary>
|
|
///
|
|
///
|
|
/// <example>
|
|
/// <code source="Unit Tests\Accord.Tests.Math\Matrix\MatrixTest.cs" region="doc_inverse" />
|
|
/// </example>
|
|
///
|
|
public static Single[][] Inverse(this Single[][] matrix)
|
|
{
|
|
return Inverse(matrix, false);
|
|
}
|
|
|
|
/// <summary>
|
|
/// Computes the inverse of a matrix.
|
|
/// </summary>
|
|
///
|
|
/// <example>
|
|
/// <code source="Unit Tests\Accord.Tests.Math\Matrix\MatrixTest.cs" region="doc_inverse" />
|
|
/// </example>
|
|
///
|
|
public static Single[][] Inverse(this Single[][] matrix, bool inPlace)
|
|
{
|
|
int rows = matrix.Length;
|
|
int cols = matrix[0].Length;
|
|
|
|
if (rows != cols)
|
|
throw new ArgumentException("Matrix must be square", "matrix");
|
|
|
|
if (rows == 3)
|
|
{
|
|
// Special case for 3x3 matrices
|
|
Single a = matrix[0][0], b = matrix[0][1], c = matrix[0][2];
|
|
Single d = matrix[1][0], e = matrix[1][1], f = matrix[1][2];
|
|
Single g = matrix[2][0], h = matrix[2][1], i = matrix[2][2];
|
|
|
|
Single den = a * (e * i - f * h) -
|
|
b * (d * i - f * g) +
|
|
c * (d * h - e * g);
|
|
|
|
if (den == 0)
|
|
throw new SingularMatrixException();
|
|
|
|
Single m = 1 / den;
|
|
|
|
var inv = matrix;
|
|
if (!inPlace)
|
|
{
|
|
inv = new Single[3][];
|
|
for (int j = 0; j < inv.Length; j++)
|
|
inv[j] = new Single[3];
|
|
}
|
|
|
|
inv[0][0] = m * (e * i - f * h);
|
|
inv[0][1] = m * (c * h - b * i);
|
|
inv[0][2] = m * (b * f - c * e);
|
|
inv[1][0] = m * (f * g - d * i);
|
|
inv[1][1] = m * (a * i - c * g);
|
|
inv[1][2] = m * (c * d - a * f);
|
|
inv[2][0] = m * (d * h - e * g);
|
|
inv[2][1] = m * (b * g - a * h);
|
|
inv[2][2] = m * (a * e - b * d);
|
|
|
|
return inv;
|
|
}
|
|
|
|
if (rows == 2)
|
|
{
|
|
// Special case for 2x2 matrices
|
|
Single a = matrix[0][0], b = matrix[0][1];
|
|
Single c = matrix[1][0], d = matrix[1][1];
|
|
|
|
Single den = a * d - b * c;
|
|
|
|
if (den == 0)
|
|
throw new SingularMatrixException();
|
|
|
|
Single m = 1 / den;
|
|
|
|
var inv = matrix;
|
|
if (!inPlace)
|
|
{
|
|
inv = new Single[2][];
|
|
for (int j = 0; j < inv.Length; j++)
|
|
inv[j] = new Single[2];
|
|
}
|
|
|
|
inv[0][0] = +m * d;
|
|
inv[0][1] = -m * b;
|
|
inv[1][0] = -m * c;
|
|
inv[1][1] = +m * a;
|
|
|
|
return inv;
|
|
}
|
|
|
|
return new JaggedLuDecompositionF(matrix, false, inPlace).Inverse();
|
|
}
|
|
|
|
/// <summary>
|
|
/// Computes the pseudo-inverse of a matrix.
|
|
/// </summary>
|
|
///
|
|
/// <example>
|
|
/// <code source="Unit Tests\Accord.Tests.Math\Matrix\MatrixTest.cs" region="doc_pseudoinverse" />
|
|
/// </example>
|
|
///
|
|
public static Single[][] PseudoInverse(this Single[][] matrix)
|
|
{
|
|
return new JaggedSingularValueDecompositionF(matrix,
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true).Inverse();
|
|
}
|
|
|
|
/// <summary>
|
|
/// Divides two matrices by multiplying A by the inverse of B.
|
|
/// </summary>
|
|
///
|
|
/// <param name="a">The first matrix.</param>
|
|
/// <param name="b">The second matrix (which will be inverted).</param>
|
|
/// <param name="leastSquares">True to produce a solution even if the
|
|
/// <paramref name="b"/> is singular; false otherwise. Default is false.</param>
|
|
///
|
|
/// <returns>The result from the division <c>AB^-1</c> of the given matrices.</returns>
|
|
///
|
|
public static Single[,] Divide(this Single[,] a, Single[,] b, bool leastSquares = false)
|
|
{
|
|
int rows = b.Rows();
|
|
int cols = b.Columns();
|
|
|
|
if (leastSquares)
|
|
{
|
|
return new SingularValueDecompositionF(b.Transpose(),
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true).Solve(a.Transpose()).Transpose();
|
|
}
|
|
|
|
|
|
if (rows == cols && cols == a.Rows())
|
|
{
|
|
// Solve by LU Decomposition if matrix is square.
|
|
return new LuDecompositionF(b, true).SolveTranspose(a);
|
|
}
|
|
else
|
|
{
|
|
if (rows <= cols)
|
|
{
|
|
// Solve by QR Decomposition if not.
|
|
return new QrDecompositionF(b, true).SolveTranspose(a);
|
|
}
|
|
else
|
|
{
|
|
return new SingularValueDecompositionF(b.Transpose(),
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true).Solve(a.Transpose()).Transpose();
|
|
}
|
|
}
|
|
}
|
|
|
|
/// <summary>
|
|
/// Divides two matrices by multiplying A by the inverse of B.
|
|
/// </summary>
|
|
///
|
|
/// <param name="a">The first matrix.</param>
|
|
/// <param name="b">The second matrix (which will be inverted).</param>
|
|
/// <param name="leastSquares">True to produce a solution even if the
|
|
/// <paramref name="b"/> is singular; false otherwise. Default is false.</param>
|
|
///
|
|
/// <returns>The result from the division <c>AB^-1</c> of the given matrices.</returns>
|
|
///
|
|
public static Single[][] Divide(this Single[][] a, Single[][] b, bool leastSquares = false)
|
|
{
|
|
int rows = b.Rows();
|
|
int cols = b.Columns();
|
|
|
|
if (leastSquares)
|
|
{
|
|
return new JaggedSingularValueDecompositionF(b.Transpose(),
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true).Solve(a.Transpose()).Transpose();
|
|
}
|
|
|
|
|
|
if (rows == cols && cols == a.Rows())
|
|
{
|
|
// Solve by LU Decomposition if matrix is square.
|
|
return new JaggedLuDecompositionF(b, true).SolveTranspose(a);
|
|
}
|
|
else
|
|
{
|
|
if (rows <= cols)
|
|
{
|
|
// Solve by QR Decomposition if not.
|
|
return new JaggedQrDecompositionF(b, true).SolveTranspose(a);
|
|
}
|
|
else
|
|
{
|
|
return new JaggedSingularValueDecompositionF(b.Transpose(),
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true).Solve(a.Transpose()).Transpose();
|
|
}
|
|
}
|
|
}
|
|
|
|
/// <summary>
|
|
/// Gets the null-space of a column vector.
|
|
/// </summary>
|
|
///
|
|
public static Single[][] Null(this Single[] vector)
|
|
{
|
|
return Null(Jagged.ColumnVector(vector));
|
|
}
|
|
|
|
/// <summary>
|
|
/// Gets the null-space of a matrix.
|
|
/// </summary>
|
|
///
|
|
public static Single[][] Null(this Single[][] matrix)
|
|
{
|
|
var qr = new JaggedQrDecompositionF(matrix, economy: false);
|
|
var Q = qr.OrthogonalFactor;
|
|
var threshold = matrix.GetLength().Max() * Constants.SingleEpsilon;
|
|
int[] idx = qr.Diagonal.Find(x => (Single)Math.Abs(x) < threshold);
|
|
return Q.GetColumns(idx);
|
|
}
|
|
|
|
/// <summary>
|
|
/// Gets the null-space of a matrix.
|
|
/// </summary>
|
|
///
|
|
public static Single[,] Null(this Single[,] matrix)
|
|
{
|
|
var qr = new QrDecompositionF(matrix, economy: false);
|
|
var Q = qr.OrthogonalFactor;
|
|
var threshold = matrix.GetLength().Max() * Constants.SingleEpsilon;
|
|
int[] idx = qr.Diagonal.Find(x => (Single)Math.Abs(x) < threshold);
|
|
return Q.GetColumns(idx);
|
|
}
|
|
|
|
/// <summary>
|
|
/// Returns the solution matrix if the matrix is square or the least squares solution otherwise.
|
|
/// </summary>
|
|
///
|
|
/// <param name="matrix">The matrix for the linear problem.</param>
|
|
/// <param name="rightSide">The right side <c>b</c>.</param>
|
|
/// <param name="leastSquares">True to produce a solution even if the
|
|
/// <paramref name="matrix"/> is singular; false otherwise. Default is false.</param>
|
|
///
|
|
/// <remarks>
|
|
/// Please note that this does not check if the matrix is non-singular
|
|
/// before attempting to solve. If a least squares solution is desired
|
|
/// in case the matrix is singular, pass true to the <paramref name="leastSquares"/>
|
|
/// parameter when calling this function.
|
|
/// </remarks>
|
|
///
|
|
/// <example>
|
|
/// <code>
|
|
/// // Create a matrix. Please note that this matrix
|
|
/// // is singular (i.e. not invertible), so only a
|
|
/// // least squares solution would be feasible here.
|
|
///
|
|
/// Decimal[,] matrix =
|
|
/// {
|
|
/// { 1, 2, 3 },
|
|
/// { 4, 5, 6 },
|
|
/// { 7, 8, 9 },
|
|
/// };
|
|
///
|
|
/// // Define a right side matrix b:
|
|
/// Decimal[,] rightSide = { {1}, {2}, {3} };
|
|
///
|
|
/// // Solve the linear system Ax = b by finding x:
|
|
/// Decimal[,] x = Matrix.Solve(matrix, rightSide, leastSquares: true);
|
|
///
|
|
/// // The answer should be { {-1/18}, {2/18}, {5/18} }.
|
|
/// </code>
|
|
/// </example>
|
|
///
|
|
public static Decimal[,] Solve(this Decimal[,] matrix, Decimal[,] rightSide, bool leastSquares = false)
|
|
{
|
|
int rows = matrix.GetLength(0);
|
|
int cols = matrix.GetLength(1);
|
|
|
|
if (rows != rightSide.GetLength(0))
|
|
throw new DimensionMismatchException("rightSide",
|
|
"The number of rows in the right hand side matrix must be "
|
|
+ "equal to the number of rows in the problem matrix.");
|
|
|
|
if (leastSquares)
|
|
{
|
|
return new SingularValueDecompositionD(matrix,
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true).Solve(rightSide);
|
|
}
|
|
|
|
|
|
if (rows == cols)
|
|
{
|
|
// Solve by LU Decomposition if matrix is square.
|
|
return new LuDecompositionD(matrix).Solve(rightSide);
|
|
}
|
|
else
|
|
{
|
|
if (cols < rows)
|
|
{
|
|
// Solve by QR Decomposition if not.
|
|
return new QrDecompositionD(matrix).Solve(rightSide);
|
|
}
|
|
else
|
|
{
|
|
return new SingularValueDecompositionD(matrix,
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true).Solve(rightSide);
|
|
}
|
|
}
|
|
}
|
|
|
|
/// <summary>
|
|
/// Returns the solution matrix if the matrix is square or the least squares solution otherwise.
|
|
/// </summary>
|
|
///
|
|
/// <param name="matrix">The matrix for the linear problem.</param>
|
|
/// <param name="rightSide">The right side <c>b</c>.</param>
|
|
/// <param name="leastSquares">True to produce a solution even if the
|
|
/// <paramref name="matrix"/> is singular; false otherwise. Default is false.</param>
|
|
///
|
|
/// <remarks>
|
|
/// Please note that this does not check if the matrix is non-singular
|
|
/// before attempting to solve. If a least squares solution is desired
|
|
/// in case the matrix is singular, pass true to the <paramref name="leastSquares"/>
|
|
/// parameter when calling this function.
|
|
/// </remarks>
|
|
///
|
|
/// <example>
|
|
/// <code>
|
|
/// // Create a matrix. Please note that this matrix
|
|
/// // is singular (i.e. not invertible), so only a
|
|
/// // least squares solution would be feasible here.
|
|
///
|
|
/// Decimal[,] matrix =
|
|
/// {
|
|
/// { 1, 2, 3 },
|
|
/// { 4, 5, 6 },
|
|
/// { 7, 8, 9 },
|
|
/// };
|
|
///
|
|
/// // Define a right side vector b:
|
|
/// Decimal[] rightSide = { 1, 2, 3 };
|
|
///
|
|
/// // Solve the linear system Ax = b by finding x:
|
|
/// Decimal[] x = Matrix.Solve(matrix, rightSide, leastSquares: true);
|
|
///
|
|
/// // The answer should be { -1/18, 2/18, 5/18 }.
|
|
/// </code>
|
|
/// </example>
|
|
///
|
|
public static Decimal[] Solve(this Decimal[,] matrix, Decimal[] rightSide, bool leastSquares = false)
|
|
{
|
|
if (matrix == null)
|
|
throw new ArgumentNullException("matrix");
|
|
|
|
if (rightSide == null)
|
|
throw new ArgumentNullException("rightSide");
|
|
|
|
int rows = matrix.GetLength(0);
|
|
int cols = matrix.GetLength(1);
|
|
|
|
if (rows != rightSide.Length)
|
|
throw new DimensionMismatchException("rightSide",
|
|
"The right hand side vector must have the same length"
|
|
+ "as there are rows of the problem matrix.");
|
|
|
|
if (leastSquares)
|
|
{
|
|
return new SingularValueDecompositionD(matrix,
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true).Solve(rightSide);
|
|
}
|
|
|
|
|
|
if (rows == cols)
|
|
{
|
|
// Solve by LU Decomposition if matrix is square.
|
|
return new LuDecompositionD(matrix).Solve(rightSide);
|
|
}
|
|
else
|
|
{
|
|
if (cols < rows)
|
|
{
|
|
// Solve by QR Decomposition if not.
|
|
return new QrDecompositionD(matrix).Solve(rightSide);
|
|
}
|
|
else
|
|
{
|
|
return new SingularValueDecompositionD(matrix,
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true).Solve(rightSide);
|
|
}
|
|
}
|
|
}
|
|
|
|
/// <summary>
|
|
/// Computes the inverse of a matrix.
|
|
/// </summary>
|
|
///
|
|
/// <example>
|
|
/// <code source="Unit Tests\Accord.Tests.Math\Matrix\MatrixTest.cs" region="doc_inverse" />
|
|
/// </example>
|
|
///
|
|
public static Decimal[,] Inverse(this Decimal[,] matrix)
|
|
{
|
|
return Inverse(matrix, false);
|
|
}
|
|
|
|
/// <summary>
|
|
/// Computes the inverse of a matrix.
|
|
/// </summary>
|
|
///
|
|
/// <example>
|
|
/// <code source="Unit Tests\Accord.Tests.Math\Matrix\MatrixTest.cs" region="doc_inverse" />
|
|
/// </example>
|
|
///
|
|
public static Decimal[,] Inverse(this Decimal[,] matrix, bool inPlace)
|
|
{
|
|
int rows = matrix.GetLength(0);
|
|
int cols = matrix.GetLength(1);
|
|
|
|
if (rows != cols)
|
|
throw new ArgumentException("Matrix must be square", "matrix");
|
|
|
|
if (rows == 3)
|
|
{
|
|
// Special case for 3x3 matrices
|
|
Decimal a = matrix[0, 0], b = matrix[0, 1], c = matrix[0, 2];
|
|
Decimal d = matrix[1, 0], e = matrix[1, 1], f = matrix[1, 2];
|
|
Decimal g = matrix[2, 0], h = matrix[2, 1], i = matrix[2, 2];
|
|
|
|
Decimal den = a * (e * i - f * h) -
|
|
b * (d * i - f * g) +
|
|
c * (d * h - e * g);
|
|
|
|
if (den == 0)
|
|
throw new SingularMatrixException();
|
|
|
|
Decimal m = 1 / den;
|
|
|
|
var inv = (inPlace) ? matrix : new Decimal[3, 3];
|
|
inv[0, 0] = m * (e * i - f * h);
|
|
inv[0, 1] = m * (c * h - b * i);
|
|
inv[0, 2] = m * (b * f - c * e);
|
|
inv[1, 0] = m * (f * g - d * i);
|
|
inv[1, 1] = m * (a * i - c * g);
|
|
inv[1, 2] = m * (c * d - a * f);
|
|
inv[2, 0] = m * (d * h - e * g);
|
|
inv[2, 1] = m * (b * g - a * h);
|
|
inv[2, 2] = m * (a * e - b * d);
|
|
|
|
return inv;
|
|
}
|
|
|
|
if (rows == 2)
|
|
{
|
|
// Special case for 2x2 matrices
|
|
Decimal a = matrix[0, 0], b = matrix[0, 1];
|
|
Decimal c = matrix[1, 0], d = matrix[1, 1];
|
|
|
|
Decimal den = a * d - b * c;
|
|
|
|
if (den == 0)
|
|
throw new SingularMatrixException();
|
|
|
|
Decimal m = 1 / den;
|
|
|
|
var inv = (inPlace) ? matrix : new Decimal[2, 2];
|
|
inv[0, 0] = +m * d;
|
|
inv[0, 1] = -m * b;
|
|
inv[1, 0] = -m * c;
|
|
inv[1, 1] = +m * a;
|
|
|
|
return inv;
|
|
}
|
|
|
|
return new LuDecompositionD(matrix, false, inPlace).Inverse();
|
|
}
|
|
|
|
/// <summary>
|
|
/// Computes the pseudo-inverse of a matrix.
|
|
/// </summary>
|
|
///
|
|
/// <example>
|
|
/// <code source="Unit Tests\Accord.Tests.Math\Matrix\MatrixTest.cs" region="doc_pseudoinverse" />
|
|
/// </example>
|
|
///
|
|
public static Decimal[,] PseudoInverse(this Decimal[,] matrix)
|
|
{
|
|
return new SingularValueDecompositionD(matrix,
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true).Inverse();
|
|
}
|
|
|
|
|
|
/// <summary>
|
|
/// Returns the solution matrix if the matrix is square or the least squares solution otherwise.
|
|
/// </summary>
|
|
///
|
|
/// <param name="matrix">The matrix for the linear problem.</param>
|
|
/// <param name="rightSide">The right side <c>b</c>.</param>
|
|
/// <param name="leastSquares">True to produce a solution even if the
|
|
/// <paramref name="matrix"/> is singular; false otherwise. Default is false.</param>
|
|
///
|
|
/// <remarks>
|
|
/// Please note that this does not check if the matrix is non-singular
|
|
/// before attempting to solve. If a least squares solution is desired
|
|
/// in case the matrix is singular, pass true to the <paramref name="leastSquares"/>
|
|
/// parameter when calling this function.
|
|
/// </remarks>
|
|
///
|
|
/// <example>
|
|
/// <code>
|
|
/// // Create a matrix. Please note that this matrix
|
|
/// // is singular (i.e. not invertible), so only a
|
|
/// // least squares solution would be feasible here.
|
|
///
|
|
/// Decimal[,] matrix =
|
|
/// {
|
|
/// { 1, 2, 3 },
|
|
/// { 4, 5, 6 },
|
|
/// { 7, 8, 9 },
|
|
/// };
|
|
///
|
|
/// // Define a right side matrix b:
|
|
/// Decimal[,] rightSide = { {1}, {2}, {3} };
|
|
///
|
|
/// // Solve the linear system Ax = b by finding x:
|
|
/// Decimal[,] x = Matrix.Solve(matrix, rightSide, leastSquares: true);
|
|
///
|
|
/// // The answer should be { {-1/18}, {2/18}, {5/18} }.
|
|
/// </code>
|
|
/// </example>
|
|
///
|
|
public static Decimal[][] Solve(this Decimal[][] matrix, Decimal[][] rightSide, bool leastSquares = false)
|
|
{
|
|
if (matrix.Length != rightSide.Length)
|
|
{
|
|
throw new DimensionMismatchException("rightSide",
|
|
"The number of rows in the right hand side matrix must be "
|
|
+ "equal to the number of rows in the problem matrix.");
|
|
}
|
|
|
|
return matrix.Decompose(leastSquares).Solve(rightSide);
|
|
}
|
|
|
|
/// <summary>
|
|
/// Returns the solution matrix if the matrix is square or the least squares solution otherwise.
|
|
/// </summary>
|
|
///
|
|
/// <param name="matrix">The matrix for the linear problem.</param>
|
|
/// <param name="rightSide">The right side <c>b</c>.</param>
|
|
/// <param name="leastSquares">True to produce a solution even if the
|
|
/// <paramref name="matrix"/> is singular; false otherwise. Default is false.</param>
|
|
///
|
|
/// <remarks>
|
|
/// Please note that this does not check if the matrix is non-singular
|
|
/// before attempting to solve. If a least squares solution is desired
|
|
/// in case the matrix is singular, pass true to the <paramref name="leastSquares"/>
|
|
/// parameter when calling this function.
|
|
/// </remarks>
|
|
///
|
|
/// <example>
|
|
/// <code>
|
|
/// // Create a matrix. Please note that this matrix
|
|
/// // is singular (i.e. not invertible), so only a
|
|
/// // least squares solution would be feasible here.
|
|
///
|
|
/// Decimal[,] matrix =
|
|
/// {
|
|
/// { 1, 2, 3 },
|
|
/// { 4, 5, 6 },
|
|
/// { 7, 8, 9 },
|
|
/// };
|
|
///
|
|
/// // Define a right side vector b:
|
|
/// Decimal[] rightSide = { 1, 2, 3 };
|
|
///
|
|
/// // Solve the linear system Ax = b by finding x:
|
|
/// Decimal[] x = Matrix.Solve(matrix, rightSide, leastSquares: true);
|
|
///
|
|
/// // The answer should be { -1/18, 2/18, 5/18 }.
|
|
/// </code>
|
|
/// </example>
|
|
///
|
|
public static Decimal[] Solve(this Decimal[][] matrix, Decimal[] rightSide, bool leastSquares = false)
|
|
{
|
|
if (matrix.Length != rightSide.Length)
|
|
{
|
|
throw new DimensionMismatchException("rightSide",
|
|
"The right hand side vector must have the same length"
|
|
+ "as there are rows of the problem matrix.");
|
|
}
|
|
|
|
return matrix.Decompose(leastSquares).Solve(rightSide);
|
|
}
|
|
|
|
/// <summary>
|
|
/// Returns the solution matrix for a linear system involving a diagonal matrix ion the right-hand side.
|
|
/// </summary>
|
|
///
|
|
/// <param name="matrix">The matrix for the linear problem.</param>
|
|
/// <param name="diagonalRightSide">The right side <c>b</c>.</param>
|
|
/// <param name="leastSquares">True to produce a solution even if the
|
|
/// <paramref name="matrix"/> is singular; false otherwise. Default is false.</param>
|
|
///
|
|
/// <remarks>
|
|
/// Please note that this does not check if the matrix is non-singular
|
|
/// before attempting to solve. If a least squares solution is desired
|
|
/// in case the matrix is singular, pass true to the <paramref name="leastSquares"/>
|
|
/// parameter when calling this function.
|
|
/// </remarks>
|
|
///
|
|
public static Decimal[][] SolveForDiagonal(this Decimal[][] matrix, Decimal[] diagonalRightSide, bool leastSquares = false)
|
|
{
|
|
if (matrix.Length != diagonalRightSide.Length)
|
|
{
|
|
throw new DimensionMismatchException("diagonalRightSide",
|
|
"The right hand side matrix must have the same length"
|
|
+ "as there are rows of the problem matrix.");
|
|
}
|
|
|
|
return matrix.Decompose(leastSquares).SolveForDiagonal(diagonalRightSide);
|
|
}
|
|
|
|
|
|
/// <summary>
|
|
/// Creates a matrix decomposition that be used to compute the solution matrix if the
|
|
/// matrix is square or the least squares solution otherwise.
|
|
/// </summary>
|
|
///
|
|
public static ISolverMatrixDecomposition<Decimal> Decompose(this Decimal[,] matrix, bool leastSquares = false)
|
|
{
|
|
int rows = matrix.Rows();
|
|
int cols = matrix.Columns();
|
|
|
|
if (leastSquares)
|
|
{
|
|
return new SingularValueDecompositionD(matrix,
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true);
|
|
}
|
|
|
|
if (rows == cols)
|
|
{
|
|
// Solve by LU Decomposition if matrix is square.
|
|
return new LuDecompositionD(matrix);
|
|
}
|
|
else
|
|
{
|
|
if (cols < rows)
|
|
{
|
|
// Solve by QR Decomposition if not.
|
|
return new QrDecompositionD(matrix);
|
|
}
|
|
else
|
|
{
|
|
return new SingularValueDecompositionD(matrix,
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true);
|
|
}
|
|
}
|
|
}
|
|
|
|
/// <summary>
|
|
/// Creates a matrix decomposition that be used to compute the solution matrix if the
|
|
/// matrix is square or the least squares solution otherwise.
|
|
/// </summary>
|
|
///
|
|
public static ISolverArrayDecomposition<Decimal> Decompose(this Decimal[][] matrix, bool leastSquares = false)
|
|
{
|
|
int rows = matrix.Rows();
|
|
int cols = matrix.Columns();
|
|
|
|
if (leastSquares)
|
|
{
|
|
return new JaggedSingularValueDecompositionD(matrix,
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true);
|
|
}
|
|
|
|
if (rows == cols)
|
|
{
|
|
// Solve by LU Decomposition if matrix is square.
|
|
return new JaggedLuDecompositionD(matrix);
|
|
}
|
|
else
|
|
{
|
|
if (cols < rows)
|
|
{
|
|
// Solve by QR Decomposition if not.
|
|
return new JaggedQrDecompositionD(matrix);
|
|
}
|
|
else
|
|
{
|
|
return new JaggedSingularValueDecompositionD(matrix,
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true);
|
|
}
|
|
}
|
|
}
|
|
|
|
/// <summary>
|
|
/// Computes the inverse of a matrix.
|
|
/// </summary>
|
|
///
|
|
///
|
|
/// <example>
|
|
/// <code source="Unit Tests\Accord.Tests.Math\Matrix\MatrixTest.cs" region="doc_inverse" />
|
|
/// </example>
|
|
///
|
|
public static Decimal[][] Inverse(this Decimal[][] matrix)
|
|
{
|
|
return Inverse(matrix, false);
|
|
}
|
|
|
|
/// <summary>
|
|
/// Computes the inverse of a matrix.
|
|
/// </summary>
|
|
///
|
|
/// <example>
|
|
/// <code source="Unit Tests\Accord.Tests.Math\Matrix\MatrixTest.cs" region="doc_inverse" />
|
|
/// </example>
|
|
///
|
|
public static Decimal[][] Inverse(this Decimal[][] matrix, bool inPlace)
|
|
{
|
|
int rows = matrix.Length;
|
|
int cols = matrix[0].Length;
|
|
|
|
if (rows != cols)
|
|
throw new ArgumentException("Matrix must be square", "matrix");
|
|
|
|
if (rows == 3)
|
|
{
|
|
// Special case for 3x3 matrices
|
|
Decimal a = matrix[0][0], b = matrix[0][1], c = matrix[0][2];
|
|
Decimal d = matrix[1][0], e = matrix[1][1], f = matrix[1][2];
|
|
Decimal g = matrix[2][0], h = matrix[2][1], i = matrix[2][2];
|
|
|
|
Decimal den = a * (e * i - f * h) -
|
|
b * (d * i - f * g) +
|
|
c * (d * h - e * g);
|
|
|
|
if (den == 0)
|
|
throw new SingularMatrixException();
|
|
|
|
Decimal m = 1 / den;
|
|
|
|
var inv = matrix;
|
|
if (!inPlace)
|
|
{
|
|
inv = new Decimal[3][];
|
|
for (int j = 0; j < inv.Length; j++)
|
|
inv[j] = new Decimal[3];
|
|
}
|
|
|
|
inv[0][0] = m * (e * i - f * h);
|
|
inv[0][1] = m * (c * h - b * i);
|
|
inv[0][2] = m * (b * f - c * e);
|
|
inv[1][0] = m * (f * g - d * i);
|
|
inv[1][1] = m * (a * i - c * g);
|
|
inv[1][2] = m * (c * d - a * f);
|
|
inv[2][0] = m * (d * h - e * g);
|
|
inv[2][1] = m * (b * g - a * h);
|
|
inv[2][2] = m * (a * e - b * d);
|
|
|
|
return inv;
|
|
}
|
|
|
|
if (rows == 2)
|
|
{
|
|
// Special case for 2x2 matrices
|
|
Decimal a = matrix[0][0], b = matrix[0][1];
|
|
Decimal c = matrix[1][0], d = matrix[1][1];
|
|
|
|
Decimal den = a * d - b * c;
|
|
|
|
if (den == 0)
|
|
throw new SingularMatrixException();
|
|
|
|
Decimal m = 1 / den;
|
|
|
|
var inv = matrix;
|
|
if (!inPlace)
|
|
{
|
|
inv = new Decimal[2][];
|
|
for (int j = 0; j < inv.Length; j++)
|
|
inv[j] = new Decimal[2];
|
|
}
|
|
|
|
inv[0][0] = +m * d;
|
|
inv[0][1] = -m * b;
|
|
inv[1][0] = -m * c;
|
|
inv[1][1] = +m * a;
|
|
|
|
return inv;
|
|
}
|
|
|
|
return new JaggedLuDecompositionD(matrix, false, inPlace).Inverse();
|
|
}
|
|
|
|
/// <summary>
|
|
/// Computes the pseudo-inverse of a matrix.
|
|
/// </summary>
|
|
///
|
|
/// <example>
|
|
/// <code source="Unit Tests\Accord.Tests.Math\Matrix\MatrixTest.cs" region="doc_pseudoinverse" />
|
|
/// </example>
|
|
///
|
|
public static Decimal[][] PseudoInverse(this Decimal[][] matrix)
|
|
{
|
|
return new JaggedSingularValueDecompositionD(matrix,
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true).Inverse();
|
|
}
|
|
|
|
/// <summary>
|
|
/// Divides two matrices by multiplying A by the inverse of B.
|
|
/// </summary>
|
|
///
|
|
/// <param name="a">The first matrix.</param>
|
|
/// <param name="b">The second matrix (which will be inverted).</param>
|
|
/// <param name="leastSquares">True to produce a solution even if the
|
|
/// <paramref name="b"/> is singular; false otherwise. Default is false.</param>
|
|
///
|
|
/// <returns>The result from the division <c>AB^-1</c> of the given matrices.</returns>
|
|
///
|
|
public static Decimal[,] Divide(this Decimal[,] a, Decimal[,] b, bool leastSquares = false)
|
|
{
|
|
int rows = b.Rows();
|
|
int cols = b.Columns();
|
|
|
|
if (leastSquares)
|
|
{
|
|
return new SingularValueDecompositionD(b.Transpose(),
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true).Solve(a.Transpose()).Transpose();
|
|
}
|
|
|
|
|
|
if (rows == cols && cols == a.Rows())
|
|
{
|
|
// Solve by LU Decomposition if matrix is square.
|
|
return new LuDecompositionD(b, true).SolveTranspose(a);
|
|
}
|
|
else
|
|
{
|
|
if (rows <= cols)
|
|
{
|
|
// Solve by QR Decomposition if not.
|
|
return new QrDecompositionD(b, true).SolveTranspose(a);
|
|
}
|
|
else
|
|
{
|
|
return new SingularValueDecompositionD(b.Transpose(),
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true).Solve(a.Transpose()).Transpose();
|
|
}
|
|
}
|
|
}
|
|
|
|
/// <summary>
|
|
/// Divides two matrices by multiplying A by the inverse of B.
|
|
/// </summary>
|
|
///
|
|
/// <param name="a">The first matrix.</param>
|
|
/// <param name="b">The second matrix (which will be inverted).</param>
|
|
/// <param name="leastSquares">True to produce a solution even if the
|
|
/// <paramref name="b"/> is singular; false otherwise. Default is false.</param>
|
|
///
|
|
/// <returns>The result from the division <c>AB^-1</c> of the given matrices.</returns>
|
|
///
|
|
public static Decimal[][] Divide(this Decimal[][] a, Decimal[][] b, bool leastSquares = false)
|
|
{
|
|
int rows = b.Rows();
|
|
int cols = b.Columns();
|
|
|
|
if (leastSquares)
|
|
{
|
|
return new JaggedSingularValueDecompositionD(b.Transpose(),
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true).Solve(a.Transpose()).Transpose();
|
|
}
|
|
|
|
|
|
if (rows == cols && cols == a.Rows())
|
|
{
|
|
// Solve by LU Decomposition if matrix is square.
|
|
return new JaggedLuDecompositionD(b, true).SolveTranspose(a);
|
|
}
|
|
else
|
|
{
|
|
if (rows <= cols)
|
|
{
|
|
// Solve by QR Decomposition if not.
|
|
return new JaggedQrDecompositionD(b, true).SolveTranspose(a);
|
|
}
|
|
else
|
|
{
|
|
return new JaggedSingularValueDecompositionD(b.Transpose(),
|
|
computeLeftSingularVectors: true,
|
|
computeRightSingularVectors: true,
|
|
autoTranspose: true).Solve(a.Transpose()).Transpose();
|
|
}
|
|
}
|
|
}
|
|
|
|
/// <summary>
|
|
/// Gets the null-space of a column vector.
|
|
/// </summary>
|
|
///
|
|
public static Decimal[][] Null(this Decimal[] vector)
|
|
{
|
|
return Null(Jagged.ColumnVector(vector));
|
|
}
|
|
|
|
/// <summary>
|
|
/// Gets the null-space of a matrix.
|
|
/// </summary>
|
|
///
|
|
public static Decimal[][] Null(this Decimal[][] matrix)
|
|
{
|
|
var qr = new JaggedQrDecompositionD(matrix, economy: false);
|
|
var Q = qr.OrthogonalFactor;
|
|
var threshold = matrix.GetLength().Max() * Constants.DecimalEpsilon;
|
|
int[] idx = qr.Diagonal.Find(x => (Decimal)Math.Abs(x) < threshold);
|
|
return Q.GetColumns(idx);
|
|
}
|
|
|
|
/// <summary>
|
|
/// Gets the null-space of a matrix.
|
|
/// </summary>
|
|
///
|
|
public static Decimal[,] Null(this Decimal[,] matrix)
|
|
{
|
|
var qr = new QrDecompositionD(matrix, economy: false);
|
|
var Q = qr.OrthogonalFactor;
|
|
var threshold = matrix.GetLength().Max() * Constants.DecimalEpsilon;
|
|
int[] idx = qr.Diagonal.Find(x => (Decimal)Math.Abs(x) < threshold);
|
|
return Q.GetColumns(idx);
|
|
}
|
|
}
|
|
} |