diff --git a/src/Numerics/LinearAlgebra/Complex/Matrix.cs b/src/Numerics/LinearAlgebra/Complex/Matrix.cs index d11ca313..21017305 100644 --- a/src/Numerics/LinearAlgebra/Complex/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Complex/Matrix.cs @@ -450,7 +450,45 @@ namespace MathNet.Numerics.LinearAlgebra.Complex Map(Complex.Log, result, Zeros.Include); } + /// + /// Computes the Moore-Penrose Pseudo-Inverse of this matrix. + /// + public override Matrix PseudoInverse() + { + var svd = Svd(true); + var w = svd.W; + var s = svd.S; + double tolerance = Math.Max(RowCount, ColumnCount) * svd.L2Norm * Precision.DoublePrecision; + for (int i = 0; i < s.Count; i++) + { + s[i] = s[i].Magnitude < tolerance ? 0 : 1/s[i]; + } + + w.SetDiagonal(s); + return (svd.U * w * svd.VT).Transpose(); + } + + /// + /// Computes the trace of this matrix. + /// + /// The trace of this matrix + /// If the matrix is not square + public override Complex Trace() + { + if (RowCount != ColumnCount) + { + throw new ArgumentException(Resources.ArgumentMatrixSquare); + } + + var sum = Complex.Zero; + for (var i = 0; i < RowCount; i++) + { + sum += At(i, i); + } + + return sum; + } /// Calculates the induced L1 norm of this matrix. /// The maximum absolute column sum of the matrix. @@ -638,27 +676,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex return Vector.Build.Dense(ret); } - /// - /// Computes the trace of this matrix. - /// - /// The trace of this matrix - /// If the matrix is not square - public override Complex Trace() - { - if (RowCount != ColumnCount) - { - throw new ArgumentException(Resources.ArgumentMatrixSquare); - } - - var sum = Complex.Zero; - for (var i = 0; i < RowCount; i++) - { - sum += At(i, i); - } - - return sum; - } - /// /// Evaluates whether this matrix is hermitian (conjugate symmetric). /// diff --git a/src/Numerics/LinearAlgebra/Complex32/Matrix.cs b/src/Numerics/LinearAlgebra/Complex32/Matrix.cs index 05c6e242..d7408613 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Matrix.cs @@ -445,6 +445,25 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 Map(Complex32.Log, result, Zeros.Include); } + /// + /// Computes the Moore-Penrose Pseudo-Inverse of this matrix. + /// + public override Matrix PseudoInverse() + { + var svd = Svd(true); + var w = svd.W; + var s = svd.S; + float tolerance = (float)(Math.Max(RowCount, ColumnCount) * svd.L2Norm * Precision.SinglePrecision); + + for (int i = 0; i < s.Count; i++) + { + s[i] = s[i].Magnitude < tolerance ? 0 : 1/s[i]; + } + + w.SetDiagonal(s); + return (svd.U * w * svd.VT).Transpose(); + } + /// Calculates the induced L1 norm of this matrix. /// The maximum absolute column sum of the matrix. public override double L1Norm() diff --git a/src/Numerics/LinearAlgebra/Double/Matrix.cs b/src/Numerics/LinearAlgebra/Double/Matrix.cs index 4f894486..4fbec2d8 100644 --- a/src/Numerics/LinearAlgebra/Double/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Double/Matrix.cs @@ -416,6 +416,46 @@ namespace MathNet.Numerics.LinearAlgebra.Double Map(Math.Log, result, Zeros.Include); } + /// + /// Computes the Moore-Penrose Pseudo-Inverse of this matrix. + /// + public override Matrix PseudoInverse() + { + var svd = Svd(true); + var w = svd.W; + var s = svd.S; + double tolerance = Math.Max(RowCount, ColumnCount) * svd.L2Norm * Precision.DoublePrecision; + + for (int i = 0; i < s.Count; i++) + { + s[i] = s[i] < tolerance ? 0 : 1/s[i]; + } + + w.SetDiagonal(s); + return (svd.U * w * svd.VT).Transpose(); + } + + /// + /// Computes the trace of this matrix. + /// + /// The trace of this matrix + /// If the matrix is not square + public override double Trace() + { + if (RowCount != ColumnCount) + { + throw new ArgumentException(Resources.ArgumentMatrixSquare); + } + + var sum = 0.0; + for (var i = 0; i < RowCount; i++) + { + sum += At(i, i); + } + + return sum; + } + /// Calculates the induced L1 norm of this matrix. /// The maximum absolute column sum of the matrix. public override double L1Norm() @@ -602,27 +642,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double return Vector.Build.Dense(ret); } - /// - /// Computes the trace of this matrix. - /// - /// The trace of this matrix - /// If the matrix is not square - public override double Trace() - { - if (RowCount != ColumnCount) - { - throw new ArgumentException(Resources.ArgumentMatrixSquare); - } - - var sum = 0.0; - for (var i = 0; i < RowCount; i++) - { - sum += At(i, i); - } - - return sum; - } - /// /// Evaluates whether this matrix is hermitian (conjugate symmetric). /// diff --git a/src/Numerics/LinearAlgebra/Matrix.Arithmetic.cs b/src/Numerics/LinearAlgebra/Matrix.Arithmetic.cs index 13ca4c02..e6b66f0d 100644 --- a/src/Numerics/LinearAlgebra/Matrix.Arithmetic.cs +++ b/src/Numerics/LinearAlgebra/Matrix.Arithmetic.cs @@ -1629,6 +1629,9 @@ namespace MathNet.Numerics.LinearAlgebra return LU().Inverse(); } + /// Computes the Moore-Penrose Pseudo-Inverse of this matrix. + public abstract Matrix PseudoInverse(); + /// /// Computes the Kronecker product of this matrix with the given matrix. The new matrix is M-by-N /// with M = this.Rows * lower.Rows and N = this.Columns * lower.Columns. diff --git a/src/Numerics/LinearAlgebra/Single/Matrix.cs b/src/Numerics/LinearAlgebra/Single/Matrix.cs index 2e1dea8b..74d86e53 100644 --- a/src/Numerics/LinearAlgebra/Single/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Single/Matrix.cs @@ -416,6 +416,25 @@ namespace MathNet.Numerics.LinearAlgebra.Single Map(x => (float)Math.Log(x), result, Zeros.Include); } + /// + /// Computes the Moore-Penrose Pseudo-Inverse of this matrix. + /// + public override Matrix PseudoInverse() + { + var svd = Svd(true); + var w = svd.W; + var s = svd.S; + float tolerance = (float)(Math.Max(RowCount, ColumnCount) * svd.L2Norm * Precision.SinglePrecision); + + for (int i = 0; i < s.Count; i++) + { + s[i] = s[i] < tolerance ? 0 : 1/s[i]; + } + + w.SetDiagonal(s); + return (svd.U * w * svd.VT).Transpose(); + } + /// /// Computes the trace of this matrix. ///