forked from tsai/mathnet-numerics
Marcus Cuda
16 years ago
97 changed files with 2719 additions and 3843 deletions
@ -0,0 +1,81 @@ |
|||
// <copyright file="Cholesky.cs" company="Math.NET">
|
|||
// Math.NET Numerics, part of the Math.NET Project
|
|||
// http://numerics.mathdotnet.com
|
|||
// http://github.com/mathnet/mathnet-numerics
|
|||
// http://mathnetnumerics.codeplex.com
|
|||
//
|
|||
// Copyright (c) 2009-2010 Math.NET
|
|||
//
|
|||
// Permission is hereby granted, free of charge, to any person
|
|||
// obtaining a copy of this software and associated documentation
|
|||
// files (the "Software"), to deal in the Software without
|
|||
// restriction, including without limitation the rights to use,
|
|||
// copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|||
// copies of the Software, and to permit persons to whom the
|
|||
// Software is furnished to do so, subject to the following
|
|||
// conditions:
|
|||
//
|
|||
// The above copyright notice and this permission notice shall be
|
|||
// included in all copies or substantial portions of the Software.
|
|||
//
|
|||
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
|||
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
|
|||
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
|||
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
|
|||
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
|
|||
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
|||
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
|
|||
// OTHER DEALINGS IN THE SOFTWARE.
|
|||
// </copyright>
|
|||
|
|||
namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization |
|||
{ |
|||
using System.Numerics; |
|||
using Generic.Factorization; |
|||
|
|||
/// <summary>
|
|||
/// <para>A class which encapsulates the functionality of a Cholesky factorization.</para>
|
|||
/// <para>For a symmetric, positive definite matrix A, the Cholesky factorization
|
|||
/// is an lower triangular matrix L so that A = L*L'.</para>
|
|||
/// </summary>
|
|||
/// <remarks>
|
|||
/// The computation of the Cholesky factorization is done at construction time. If the matrix is not symmetric
|
|||
/// or positive definite, the constructor will throw an exception.
|
|||
/// </remarks>
|
|||
public abstract class Cholesky : Cholesky<Complex> |
|||
{ |
|||
/// <summary>
|
|||
/// Gets the determinant of the matrix for which the Cholesky matrix was computed.
|
|||
/// </summary>
|
|||
public override Complex Determinant |
|||
{ |
|||
get |
|||
{ |
|||
var det = Complex.One; |
|||
for (var j = 0; j < CholeskyFactor.RowCount; j++) |
|||
{ |
|||
det *= CholeskyFactor[j, j] * CholeskyFactor[j, j]; |
|||
} |
|||
|
|||
return det; |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Gets the log determinant of the matrix for which the Cholesky matrix was computed.
|
|||
/// </summary>
|
|||
public override Complex DeterminantLn |
|||
{ |
|||
get |
|||
{ |
|||
var det = Complex.Zero; |
|||
for (var j = 0; j < CholeskyFactor.RowCount; j++) |
|||
{ |
|||
det += 2.0 * CholeskyFactor[j, j].NaturalLogarithm(); |
|||
} |
|||
|
|||
return det; |
|||
} |
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,114 @@ |
|||
// <copyright file="Evd.cs" company="Math.NET">
|
|||
// Math.NET Numerics, part of the Math.NET Project
|
|||
// http://numerics.mathdotnet.com
|
|||
// http://github.com/mathnet/mathnet-numerics
|
|||
// http://mathnetnumerics.codeplex.com
|
|||
// Copyright (c) 2009-2010 Math.NET
|
|||
// Permission is hereby granted, free of charge, to any person
|
|||
// obtaining a copy of this software and associated documentation
|
|||
// files (the "Software"), to deal in the Software without
|
|||
// restriction, including without limitation the rights to use,
|
|||
// copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|||
// copies of the Software, and to permit persons to whom the
|
|||
// Software is furnished to do so, subject to the following
|
|||
// conditions:
|
|||
// The above copyright notice and this permission notice shall be
|
|||
// included in all copies or substantial portions of the Software.
|
|||
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
|||
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
|
|||
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
|||
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
|
|||
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
|
|||
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
|||
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
|
|||
// OTHER DEALINGS IN THE SOFTWARE.
|
|||
// </copyright>
|
|||
|
|||
namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization |
|||
{ |
|||
using System.Numerics; |
|||
using Generic.Factorization; |
|||
|
|||
/// <summary>
|
|||
/// Eigenvalues and eigenvectors of a real matrix.
|
|||
/// </summary>
|
|||
/// <remarks>
|
|||
/// If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is
|
|||
/// diagonal and the eigenvector matrix V is orthogonal.
|
|||
/// I.e. A = V*D*V' and V*VT=I.
|
|||
/// If A is not symmetric, then the eigenvalue matrix D is block diagonal
|
|||
/// with the real eigenvalues in 1-by-1 blocks and any complex eigenvalues,
|
|||
/// lambda + i*mu, in 2-by-2 blocks, [lambda, mu; -mu, lambda]. The
|
|||
/// columns of V represent the eigenvectors in the sense that A*V = V*D,
|
|||
/// i.e. A.Multiply(V) equals V.Multiply(D). The matrix V may be badly
|
|||
/// conditioned, or even singular, so the validity of the equation
|
|||
/// A = V*D*Inverse(V) depends upon V.Condition().
|
|||
/// </remarks>
|
|||
public abstract class Evd : Evd<Complex> |
|||
{ |
|||
/// <summary>
|
|||
/// Gets the absolute value of determinant of the square matrix for which the EVD was computed.
|
|||
/// </summary>
|
|||
public override Complex Determinant |
|||
{ |
|||
get |
|||
{ |
|||
var det = Complex.One; |
|||
for (var i = 0; i < VectorEv.Count; i++) |
|||
{ |
|||
det *= VectorEv[i]; |
|||
|
|||
if (VectorEv[i].AlmostEqual(Complex.Zero)) |
|||
{ |
|||
return 0; |
|||
} |
|||
} |
|||
|
|||
return det.Magnitude; |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Gets the effective numerical matrix rank.
|
|||
/// </summary>
|
|||
/// <value>The number of non-negligible singular values.</value>
|
|||
public override int Rank |
|||
{ |
|||
get |
|||
{ |
|||
var rank = 0; |
|||
for (var i = 0; i < VectorEv.Count; i++) |
|||
{ |
|||
if (VectorEv[i].AlmostEqual(Complex.Zero)) |
|||
{ |
|||
continue; |
|||
} |
|||
|
|||
rank++; |
|||
} |
|||
|
|||
return rank; |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Gets a value indicating whether the matrix is full rank or not.
|
|||
/// </summary>
|
|||
/// <value><c>true</c> if the matrix is full rank; otherwise <c>false</c>.</value>
|
|||
public override bool IsFullRank |
|||
{ |
|||
get |
|||
{ |
|||
for (var i = 0; i < VectorEv.Count; i++) |
|||
{ |
|||
if (VectorEv[i].AlmostEqual(Complex.Zero)) |
|||
{ |
|||
return false; |
|||
} |
|||
} |
|||
|
|||
return true; |
|||
} |
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,89 @@ |
|||
// <copyright file="GramSchmidt.cs" company="Math.NET">
|
|||
// Math.NET Numerics, part of the Math.NET Project
|
|||
// http://numerics.mathdotnet.com
|
|||
// http://github.com/mathnet/mathnet-numerics
|
|||
// http://mathnetnumerics.codeplex.com
|
|||
// Copyright (c) 2009-2010 Math.NET
|
|||
// Permission is hereby granted, free of charge, to any person
|
|||
// obtaining a copy of this software and associated documentation
|
|||
// files (the "Software"), to deal in the Software without
|
|||
// restriction, including without limitation the rights to use,
|
|||
// copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|||
// copies of the Software, and to permit persons to whom the
|
|||
// Software is furnished to do so, subject to the following
|
|||
// conditions:
|
|||
// The above copyright notice and this permission notice shall be
|
|||
// included in all copies or substantial portions of the Software.
|
|||
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
|||
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
|
|||
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
|||
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
|
|||
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
|
|||
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
|||
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
|
|||
// OTHER DEALINGS IN THE SOFTWARE.
|
|||
// </copyright>
|
|||
|
|||
namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization |
|||
{ |
|||
using System; |
|||
using System.Numerics; |
|||
using Generic.Factorization; |
|||
using Properties; |
|||
|
|||
/// <summary>
|
|||
/// <para>A class which encapsulates the functionality of the QR decomposition Modified Gram-Schmidt Orthogonalization.</para>
|
|||
/// <para>Any real square matrix A may be decomposed as A = QR where Q is an orthogonal mxn matrix and R is an nxn upper triangular matrix.</para>
|
|||
/// </summary>
|
|||
/// <remarks>
|
|||
/// The computation of the QR decomposition is done at construction time by modified Gram-Schmidt Orthogonalization.
|
|||
/// </remarks>
|
|||
public abstract class GramSchmidt : GramSchmidt<Complex> |
|||
{ |
|||
/// <summary>
|
|||
/// Gets the absolute determinant value of the matrix for which the QR matrix was computed.
|
|||
/// </summary>
|
|||
public override Complex Determinant |
|||
{ |
|||
get |
|||
{ |
|||
if (MatrixR.RowCount != MatrixR.ColumnCount) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentMatrixSquare); |
|||
} |
|||
|
|||
var det = Complex.One; |
|||
for (var i = 0; i < MatrixR.ColumnCount; i++) |
|||
{ |
|||
det *= MatrixR.At(i, i); |
|||
if (MatrixR.At(i, i).Magnitude.AlmostEqual(0.0)) |
|||
{ |
|||
return 0; |
|||
} |
|||
} |
|||
|
|||
return det.Magnitude; |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Gets a value indicating whether the matrix is full rank or not.
|
|||
/// </summary>
|
|||
/// <value><c>true</c> if the matrix is full rank; otherwise <c>false</c>.</value>
|
|||
public override bool IsFullRank |
|||
{ |
|||
get |
|||
{ |
|||
for (var i = 0; i < MatrixR.ColumnCount; i++) |
|||
{ |
|||
if (MatrixR.At(i, i).Magnitude.AlmostEqual(0.0)) |
|||
{ |
|||
return false; |
|||
} |
|||
} |
|||
|
|||
return true; |
|||
} |
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,68 @@ |
|||
// <copyright file="LU.cs" company="Math.NET">
|
|||
// Math.NET Numerics, part of the Math.NET Project
|
|||
// http://numerics.mathdotnet.com
|
|||
// http://github.com/mathnet/mathnet-numerics
|
|||
// http://mathnetnumerics.codeplex.com
|
|||
// Copyright (c) 2009-2010 Math.NET
|
|||
// Permission is hereby granted, free of charge, to any person
|
|||
// obtaining a copy of this software and associated documentation
|
|||
// files (the "Software"), to deal in the Software without
|
|||
// restriction, including without limitation the rights to use,
|
|||
// copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|||
// copies of the Software, and to permit persons to whom the
|
|||
// Software is furnished to do so, subject to the following
|
|||
// conditions:
|
|||
// The above copyright notice and this permission notice shall be
|
|||
// included in all copies or substantial portions of the Software.
|
|||
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
|||
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
|
|||
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
|||
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
|
|||
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
|
|||
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
|||
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
|
|||
// OTHER DEALINGS IN THE SOFTWARE.
|
|||
// </copyright>
|
|||
|
|||
namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization |
|||
{ |
|||
using System.Numerics; |
|||
using Generic.Factorization; |
|||
|
|||
/// <summary>
|
|||
/// <para>A class which encapsulates the functionality of an LU factorization.</para>
|
|||
/// <para>For a matrix A, the LU factorization is a pair of lower triangular matrix L and
|
|||
/// upper triangular matrix U so that A = L*U.</para>
|
|||
/// <para>In the Math.Net implementation we also store a set of pivot elements for increased
|
|||
/// numerical stability. The pivot elements encode a permutation matrix P such that P*A = L*U.</para>
|
|||
/// </summary>
|
|||
/// <remarks>
|
|||
/// The computation of the LU factorization is done at construction time.
|
|||
/// </remarks>
|
|||
public abstract class LU : LU<Complex> |
|||
{ |
|||
/// <summary>
|
|||
/// Gets the determinant of the matrix for which the LU factorization was computed.
|
|||
/// </summary>
|
|||
public override Complex Determinant |
|||
{ |
|||
get |
|||
{ |
|||
var det = Complex.One; |
|||
for (var j = 0; j < Factors.RowCount; j++) |
|||
{ |
|||
if (Pivots[j] != j) |
|||
{ |
|||
det *= -Factors.At(j, j); |
|||
} |
|||
else |
|||
{ |
|||
det *= Factors.At(j, j); |
|||
} |
|||
} |
|||
|
|||
return det; |
|||
} |
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,91 @@ |
|||
// <copyright file="QR.cs" company="Math.NET">
|
|||
// Math.NET Numerics, part of the Math.NET Project
|
|||
// http://numerics.mathdotnet.com
|
|||
// http://github.com/mathnet/mathnet-numerics
|
|||
// http://mathnetnumerics.codeplex.com
|
|||
// Copyright (c) 2009-2010 Math.NET
|
|||
// Permission is hereby granted, free of charge, to any person
|
|||
// obtaining a copy of this software and associated documentation
|
|||
// files (the "Software"), to deal in the Software without
|
|||
// restriction, including without limitation the rights to use,
|
|||
// copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|||
// copies of the Software, and to permit persons to whom the
|
|||
// Software is furnished to do so, subject to the following
|
|||
// conditions:
|
|||
// The above copyright notice and this permission notice shall be
|
|||
// included in all copies or substantial portions of the Software.
|
|||
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
|||
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
|
|||
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
|||
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
|
|||
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
|
|||
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
|||
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
|
|||
// OTHER DEALINGS IN THE SOFTWARE.
|
|||
// </copyright>
|
|||
|
|||
namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization |
|||
{ |
|||
using System; |
|||
using System.Numerics; |
|||
using Generic.Factorization; |
|||
using Properties; |
|||
|
|||
/// <summary>
|
|||
/// <para>A class which encapsulates the functionality of the QR decomposition.</para>
|
|||
/// <para>Any real square matrix A (m x n) may be decomposed as A = QR where Q is an orthogonal matrix (m x m)
|
|||
/// (its columns are orthogonal unit vectors meaning QTQ = I) and R (m x n) is an upper triangular matrix
|
|||
/// (also called right triangular matrix).</para>
|
|||
/// </summary>
|
|||
/// <remarks>
|
|||
/// The computation of the QR decomposition is done at construction time by Householder transformation.
|
|||
/// </remarks>
|
|||
public abstract class QR : QR<Complex> |
|||
{ |
|||
/// <summary>
|
|||
/// Gets the absolute determinant value of the matrix for which the QR matrix was computed.
|
|||
/// </summary>
|
|||
public override Complex Determinant |
|||
{ |
|||
get |
|||
{ |
|||
if (MatrixR.RowCount != MatrixR.ColumnCount) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentMatrixSquare); |
|||
} |
|||
|
|||
var det = Complex.One; |
|||
for (var i = 0; i < MatrixR.ColumnCount; i++) |
|||
{ |
|||
det *= MatrixR.At(i, i); |
|||
if (MatrixR.At(i, i).Magnitude.AlmostEqual(0.0)) |
|||
{ |
|||
return 0; |
|||
} |
|||
} |
|||
|
|||
return det.Magnitude; |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Gets a value indicating whether the matrix is full rank or not.
|
|||
/// </summary>
|
|||
/// <value><c>true</c> if the matrix is full rank; otherwise <c>false</c>.</value>
|
|||
public override bool IsFullRank |
|||
{ |
|||
get |
|||
{ |
|||
for (var i = 0; i < MatrixR.ColumnCount; i++) |
|||
{ |
|||
if (MatrixR.At(i, i).Magnitude.AlmostEqual(0.0)) |
|||
{ |
|||
return false; |
|||
} |
|||
} |
|||
|
|||
return true; |
|||
} |
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,119 @@ |
|||
// <copyright file="Svd.cs" company="Math.NET">
|
|||
// Math.NET Numerics, part of the Math.NET Project
|
|||
// http://numerics.mathdotnet.com
|
|||
// http://github.com/mathnet/mathnet-numerics
|
|||
// http://mathnetnumerics.codeplex.com
|
|||
//
|
|||
// Copyright (c) 2009-2010 Math.NET
|
|||
//
|
|||
// Permission is hereby granted, free of charge, to any person
|
|||
// obtaining a copy of this software and associated documentation
|
|||
// files (the "Software"), to deal in the Software without
|
|||
// restriction, including without limitation the rights to use,
|
|||
// copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|||
// copies of the Software, and to permit persons to whom the
|
|||
// Software is furnished to do so, subject to the following
|
|||
// conditions:
|
|||
//
|
|||
// The above copyright notice and this permission notice shall be
|
|||
// included in all copies or substantial portions of the Software.
|
|||
//
|
|||
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
|||
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
|
|||
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
|||
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
|
|||
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
|
|||
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
|||
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
|
|||
// OTHER DEALINGS IN THE SOFTWARE.
|
|||
// </copyright>
|
|||
|
|||
namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization |
|||
{ |
|||
using System; |
|||
using System.Linq; |
|||
using System.Numerics; |
|||
using Generic; |
|||
using Generic.Factorization; |
|||
using Properties; |
|||
|
|||
/// <summary>
|
|||
/// <para>A class which encapsulates the functionality of the singular value decomposition (SVD).</para>
|
|||
/// <para>Suppose M is an m-by-n matrix whose entries are real numbers.
|
|||
/// Then there exists a factorization of the form M = UΣVT where:
|
|||
/// - U is an m-by-m unitary matrix;
|
|||
/// - Σ is m-by-n diagonal matrix with nonnegative real numbers on the diagonal;
|
|||
/// - VT denotes transpose of V, an n-by-n unitary matrix;
|
|||
/// Such a factorization is called a singular-value decomposition of M. A common convention is to order the diagonal
|
|||
/// entries Σ(i,i) in descending order. In this case, the diagonal matrix Σ is uniquely determined
|
|||
/// by M (though the matrices U and V are not). The diagonal entries of Σ are known as the singular values of M.</para>
|
|||
/// </summary>
|
|||
/// <remarks>
|
|||
/// The computation of the singular value decomposition is done at construction time.
|
|||
/// </remarks>
|
|||
public abstract class Svd : Svd<Complex> |
|||
{ |
|||
/// <summary>
|
|||
/// Gets the effective numerical matrix rank.
|
|||
/// </summary>
|
|||
/// <value>The number of non-negligible singular values.</value>
|
|||
public override int Rank |
|||
{ |
|||
get |
|||
{ |
|||
return VectorS.Count(t => !t.Magnitude.AlmostEqual(0.0)); |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Gets the two norm of the <see cref="Matrix{T}"/>.
|
|||
/// </summary>
|
|||
/// <returns>The 2-norm of the <see cref="Matrix{T}"/>.</returns>
|
|||
public override Complex Norm2 |
|||
{ |
|||
get |
|||
{ |
|||
return VectorS[0].Magnitude; |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Gets the condition number <b>max(S) / min(S)</b>
|
|||
/// </summary>
|
|||
/// <returns>The condition number.</returns>
|
|||
public override Complex ConditionNumber |
|||
{ |
|||
get |
|||
{ |
|||
var tmp = Math.Min(MatrixU.RowCount, MatrixVT.ColumnCount) - 1; |
|||
return VectorS[0].Magnitude / VectorS[tmp].Magnitude; |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Gets the determinant of the square matrix for which the SVD was computed.
|
|||
/// </summary>
|
|||
public override Complex Determinant |
|||
{ |
|||
get |
|||
{ |
|||
if (MatrixU.RowCount != MatrixVT.ColumnCount) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentMatrixSquare); |
|||
} |
|||
|
|||
var det = Complex.One; |
|||
foreach (var value in VectorS) |
|||
{ |
|||
det *= value; |
|||
if (value.Magnitude.AlmostEqual(0.0)) |
|||
{ |
|||
return 0; |
|||
} |
|||
} |
|||
|
|||
return det.Magnitude; |
|||
} |
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,81 @@ |
|||
// <copyright file="Cholesky.cs" company="Math.NET">
|
|||
// Math.NET Numerics, part of the Math.NET Project
|
|||
// http://numerics.mathdotnet.com
|
|||
// http://github.com/mathnet/mathnet-numerics
|
|||
// http://mathnetnumerics.codeplex.com
|
|||
//
|
|||
// Copyright (c) 2009-2010 Math.NET
|
|||
//
|
|||
// Permission is hereby granted, free of charge, to any person
|
|||
// obtaining a copy of this software and associated documentation
|
|||
// files (the "Software"), to deal in the Software without
|
|||
// restriction, including without limitation the rights to use,
|
|||
// copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|||
// copies of the Software, and to permit persons to whom the
|
|||
// Software is furnished to do so, subject to the following
|
|||
// conditions:
|
|||
//
|
|||
// The above copyright notice and this permission notice shall be
|
|||
// included in all copies or substantial portions of the Software.
|
|||
//
|
|||
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
|||
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
|
|||
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
|||
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
|
|||
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
|
|||
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
|||
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
|
|||
// OTHER DEALINGS IN THE SOFTWARE.
|
|||
// </copyright>
|
|||
|
|||
namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|||
{ |
|||
using Generic.Factorization; |
|||
using Numerics; |
|||
|
|||
/// <summary>
|
|||
/// <para>A class which encapsulates the functionality of a Cholesky factorization.</para>
|
|||
/// <para>For a symmetric, positive definite matrix A, the Cholesky factorization
|
|||
/// is an lower triangular matrix L so that A = L*L'.</para>
|
|||
/// </summary>
|
|||
/// <remarks>
|
|||
/// The computation of the Cholesky factorization is done at construction time. If the matrix is not symmetric
|
|||
/// or positive definite, the constructor will throw an exception.
|
|||
/// </remarks>
|
|||
public abstract class Cholesky : Cholesky<Complex32> |
|||
{ |
|||
/// <summary>
|
|||
/// Gets the determinant of the matrix for which the Cholesky matrix was computed.
|
|||
/// </summary>
|
|||
public override Complex32 Determinant |
|||
{ |
|||
get |
|||
{ |
|||
var det = Complex32.One; |
|||
for (var j = 0; j < CholeskyFactor.RowCount; j++) |
|||
{ |
|||
det *= CholeskyFactor[j, j] * CholeskyFactor[j, j]; |
|||
} |
|||
|
|||
return det; |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Gets the log determinant of the matrix for which the Cholesky matrix was computed.
|
|||
/// </summary>
|
|||
public override Complex32 DeterminantLn |
|||
{ |
|||
get |
|||
{ |
|||
var det = Complex32.Zero; |
|||
for (var j = 0; j < CholeskyFactor.RowCount; j++) |
|||
{ |
|||
det += 2.0f * CholeskyFactor[j, j].NaturalLogarithm(); |
|||
} |
|||
|
|||
return det; |
|||
} |
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,116 @@ |
|||
// <copyright file="Evd.cs" company="Math.NET">
|
|||
// Math.NET Numerics, part of the Math.NET Project
|
|||
// http://numerics.mathdotnet.com
|
|||
// http://github.com/mathnet/mathnet-numerics
|
|||
// http://mathnetnumerics.codeplex.com
|
|||
// Copyright (c) 2009-2010 Math.NET
|
|||
// Permission is hereby granted, free of charge, to any person
|
|||
// obtaining a copy of this software and associated documentation
|
|||
// files (the "Software"), to deal in the Software without
|
|||
// restriction, including without limitation the rights to use,
|
|||
// copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|||
// copies of the Software, and to permit persons to whom the
|
|||
// Software is furnished to do so, subject to the following
|
|||
// conditions:
|
|||
// The above copyright notice and this permission notice shall be
|
|||
// included in all copies or substantial portions of the Software.
|
|||
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
|||
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
|
|||
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
|||
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
|
|||
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
|
|||
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
|||
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
|
|||
// OTHER DEALINGS IN THE SOFTWARE.
|
|||
// </copyright>
|
|||
|
|||
namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|||
{ |
|||
using System; |
|||
using System.Numerics; |
|||
using Generic.Factorization; |
|||
using Numerics; |
|||
|
|||
/// <summary>
|
|||
/// Eigenvalues and eigenvectors of a real matrix.
|
|||
/// </summary>
|
|||
/// <remarks>
|
|||
/// If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is
|
|||
/// diagonal and the eigenvector matrix V is orthogonal.
|
|||
/// I.e. A = V*D*V' and V*VT=I.
|
|||
/// If A is not symmetric, then the eigenvalue matrix D is block diagonal
|
|||
/// with the real eigenvalues in 1-by-1 blocks and any complex eigenvalues,
|
|||
/// lambda + i*mu, in 2-by-2 blocks, [lambda, mu; -mu, lambda]. The
|
|||
/// columns of V represent the eigenvectors in the sense that A*V = V*D,
|
|||
/// i.e. A.Multiply(V) equals V.Multiply(D). The matrix V may be badly
|
|||
/// conditioned, or even singular, so the validity of the equation
|
|||
/// A = V*D*Inverse(V) depends upon V.Condition().
|
|||
/// </remarks>
|
|||
public abstract class Evd : Evd<Complex32> |
|||
{ |
|||
/// <summary>
|
|||
/// Gets the absolute value of determinant of the square matrix for which the EVD was computed.
|
|||
/// </summary>
|
|||
public override Complex32 Determinant |
|||
{ |
|||
get |
|||
{ |
|||
var det = Complex.One; |
|||
for (var i = 0; i < VectorEv.Count; i++) |
|||
{ |
|||
det *= VectorEv[i]; |
|||
|
|||
if (((Complex32)VectorEv[i]).AlmostEqual(Complex32.Zero)) |
|||
{ |
|||
return 0; |
|||
} |
|||
} |
|||
|
|||
return new Complex32(Convert.ToSingle(det.Magnitude), 0.0f); |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Gets the effective numerical matrix rank.
|
|||
/// </summary>
|
|||
/// <value>The number of non-negligible singular values.</value>
|
|||
public override int Rank |
|||
{ |
|||
get |
|||
{ |
|||
var rank = 0; |
|||
for (var i = 0; i < VectorEv.Count; i++) |
|||
{ |
|||
if (((Complex32)VectorEv[i]).AlmostEqual(Complex32.Zero)) |
|||
{ |
|||
continue; |
|||
} |
|||
|
|||
rank++; |
|||
} |
|||
|
|||
return rank; |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Gets a value indicating whether the matrix is full rank or not.
|
|||
/// </summary>
|
|||
/// <value><c>true</c> if the matrix is full rank; otherwise <c>false</c>.</value>
|
|||
public override bool IsFullRank |
|||
{ |
|||
get |
|||
{ |
|||
for (var i = 0; i < VectorEv.Count; i++) |
|||
{ |
|||
if (VectorEv[i].AlmostEqual(Complex.Zero)) |
|||
{ |
|||
return false; |
|||
} |
|||
} |
|||
|
|||
return true; |
|||
} |
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,89 @@ |
|||
// <copyright file="GramSchmidt.cs" company="Math.NET">
|
|||
// Math.NET Numerics, part of the Math.NET Project
|
|||
// http://numerics.mathdotnet.com
|
|||
// http://github.com/mathnet/mathnet-numerics
|
|||
// http://mathnetnumerics.codeplex.com
|
|||
// Copyright (c) 2009-2010 Math.NET
|
|||
// Permission is hereby granted, free of charge, to any person
|
|||
// obtaining a copy of this software and associated documentation
|
|||
// files (the "Software"), to deal in the Software without
|
|||
// restriction, including without limitation the rights to use,
|
|||
// copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|||
// copies of the Software, and to permit persons to whom the
|
|||
// Software is furnished to do so, subject to the following
|
|||
// conditions:
|
|||
// The above copyright notice and this permission notice shall be
|
|||
// included in all copies or substantial portions of the Software.
|
|||
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
|||
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
|
|||
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
|||
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
|
|||
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
|
|||
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
|||
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
|
|||
// OTHER DEALINGS IN THE SOFTWARE.
|
|||
// </copyright>
|
|||
|
|||
namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|||
{ |
|||
using System; |
|||
using Generic.Factorization; |
|||
using Numerics; |
|||
using Properties; |
|||
|
|||
/// <summary>
|
|||
/// <para>A class which encapsulates the functionality of the QR decomposition Modified Gram-Schmidt Orthogonalization.</para>
|
|||
/// <para>Any real square matrix A may be decomposed as A = QR where Q is an orthogonal mxn matrix and R is an nxn upper triangular matrix.</para>
|
|||
/// </summary>
|
|||
/// <remarks>
|
|||
/// The computation of the QR decomposition is done at construction time by modified Gram-Schmidt Orthogonalization.
|
|||
/// </remarks>
|
|||
public abstract class GramSchmidt : GramSchmidt<Complex32> |
|||
{ |
|||
/// <summary>
|
|||
/// Gets the absolute determinant value of the matrix for which the QR matrix was computed.
|
|||
/// </summary>
|
|||
public override Complex32 Determinant |
|||
{ |
|||
get |
|||
{ |
|||
if (MatrixR.RowCount != MatrixR.ColumnCount) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentMatrixSquare); |
|||
} |
|||
|
|||
var det = Complex32.One; |
|||
for (var i = 0; i < MatrixR.ColumnCount; i++) |
|||
{ |
|||
det *= MatrixR.At(i, i); |
|||
if (MatrixR.At(i, i).Magnitude.AlmostEqual(0.0f)) |
|||
{ |
|||
return 0; |
|||
} |
|||
} |
|||
|
|||
return det.Magnitude; |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Gets a value indicating whether the matrix is full rank or not.
|
|||
/// </summary>
|
|||
/// <value><c>true</c> if the matrix is full rank; otherwise <c>false</c>.</value>
|
|||
public override bool IsFullRank |
|||
{ |
|||
get |
|||
{ |
|||
for (var i = 0; i < MatrixR.ColumnCount; i++) |
|||
{ |
|||
if (MatrixR.At(i, i).Magnitude.AlmostEqual(0.0f)) |
|||
{ |
|||
return false; |
|||
} |
|||
} |
|||
|
|||
return true; |
|||
} |
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,68 @@ |
|||
// <copyright file="LU.cs" company="Math.NET">
|
|||
// Math.NET Numerics, part of the Math.NET Project
|
|||
// http://numerics.mathdotnet.com
|
|||
// http://github.com/mathnet/mathnet-numerics
|
|||
// http://mathnetnumerics.codeplex.com
|
|||
// Copyright (c) 2009-2010 Math.NET
|
|||
// Permission is hereby granted, free of charge, to any person
|
|||
// obtaining a copy of this software and associated documentation
|
|||
// files (the "Software"), to deal in the Software without
|
|||
// restriction, including without limitation the rights to use,
|
|||
// copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|||
// copies of the Software, and to permit persons to whom the
|
|||
// Software is furnished to do so, subject to the following
|
|||
// conditions:
|
|||
// The above copyright notice and this permission notice shall be
|
|||
// included in all copies or substantial portions of the Software.
|
|||
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
|||
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
|
|||
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
|||
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
|
|||
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
|
|||
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
|||
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
|
|||
// OTHER DEALINGS IN THE SOFTWARE.
|
|||
// </copyright>
|
|||
|
|||
namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|||
{ |
|||
using Generic.Factorization; |
|||
using Numerics; |
|||
|
|||
/// <summary>
|
|||
/// <para>A class which encapsulates the functionality of an LU factorization.</para>
|
|||
/// <para>For a matrix A, the LU factorization is a pair of lower triangular matrix L and
|
|||
/// upper triangular matrix U so that A = L*U.</para>
|
|||
/// <para>In the Math.Net implementation we also store a set of pivot elements for increased
|
|||
/// numerical stability. The pivot elements encode a permutation matrix P such that P*A = L*U.</para>
|
|||
/// </summary>
|
|||
/// <remarks>
|
|||
/// The computation of the LU factorization is done at construction time.
|
|||
/// </remarks>
|
|||
public abstract class LU : LU<Complex32> |
|||
{ |
|||
/// <summary>
|
|||
/// Gets the determinant of the matrix for which the LU factorization was computed.
|
|||
/// </summary>
|
|||
public override Complex32 Determinant |
|||
{ |
|||
get |
|||
{ |
|||
var det = Complex32.One; |
|||
for (var j = 0; j < Factors.RowCount; j++) |
|||
{ |
|||
if (Pivots[j] != j) |
|||
{ |
|||
det *= -Factors.At(j, j); |
|||
} |
|||
else |
|||
{ |
|||
det *= Factors.At(j, j); |
|||
} |
|||
} |
|||
|
|||
return det; |
|||
} |
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,91 @@ |
|||
// <copyright file="QR.cs" company="Math.NET">
|
|||
// Math.NET Numerics, part of the Math.NET Project
|
|||
// http://numerics.mathdotnet.com
|
|||
// http://github.com/mathnet/mathnet-numerics
|
|||
// http://mathnetnumerics.codeplex.com
|
|||
// Copyright (c) 2009-2010 Math.NET
|
|||
// Permission is hereby granted, free of charge, to any person
|
|||
// obtaining a copy of this software and associated documentation
|
|||
// files (the "Software"), to deal in the Software without
|
|||
// restriction, including without limitation the rights to use,
|
|||
// copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|||
// copies of the Software, and to permit persons to whom the
|
|||
// Software is furnished to do so, subject to the following
|
|||
// conditions:
|
|||
// The above copyright notice and this permission notice shall be
|
|||
// included in all copies or substantial portions of the Software.
|
|||
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
|||
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
|
|||
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
|||
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
|
|||
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
|
|||
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
|||
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
|
|||
// OTHER DEALINGS IN THE SOFTWARE.
|
|||
// </copyright>
|
|||
|
|||
namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|||
{ |
|||
using System; |
|||
using Generic.Factorization; |
|||
using Numerics; |
|||
using Properties; |
|||
|
|||
/// <summary>
|
|||
/// <para>A class which encapsulates the functionality of the QR decomposition.</para>
|
|||
/// <para>Any real square matrix A (m x n) may be decomposed as A = QR where Q is an orthogonal matrix (m x m)
|
|||
/// (its columns are orthogonal unit vectors meaning QTQ = I) and R (m x n) is an upper triangular matrix
|
|||
/// (also called right triangular matrix).</para>
|
|||
/// </summary>
|
|||
/// <remarks>
|
|||
/// The computation of the QR decomposition is done at construction time by Householder transformation.
|
|||
/// </remarks>
|
|||
public abstract class QR : QR<Complex32> |
|||
{ |
|||
/// <summary>
|
|||
/// Gets the absolute determinant value of the matrix for which the QR matrix was computed.
|
|||
/// </summary>
|
|||
public override Complex32 Determinant |
|||
{ |
|||
get |
|||
{ |
|||
if (MatrixR.RowCount != MatrixR.ColumnCount) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentMatrixSquare); |
|||
} |
|||
|
|||
var det = Complex32.One; |
|||
for (var i = 0; i < MatrixR.ColumnCount; i++) |
|||
{ |
|||
det *= MatrixR.At(i, i); |
|||
if (MatrixR.At(i, i).Magnitude.AlmostEqual(0.0f)) |
|||
{ |
|||
return 0; |
|||
} |
|||
} |
|||
|
|||
return det.Magnitude; |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Gets a value indicating whether the matrix is full rank or not.
|
|||
/// </summary>
|
|||
/// <value><c>true</c> if the matrix is full rank; otherwise <c>false</c>.</value>
|
|||
public override bool IsFullRank |
|||
{ |
|||
get |
|||
{ |
|||
for (var i = 0; i < MatrixR.ColumnCount; i++) |
|||
{ |
|||
if (MatrixR.At(i, i).Magnitude.AlmostEqual(0.0f)) |
|||
{ |
|||
return false; |
|||
} |
|||
} |
|||
|
|||
return true; |
|||
} |
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,119 @@ |
|||
// <copyright file="Svd.cs" company="Math.NET">
|
|||
// Math.NET Numerics, part of the Math.NET Project
|
|||
// http://numerics.mathdotnet.com
|
|||
// http://github.com/mathnet/mathnet-numerics
|
|||
// http://mathnetnumerics.codeplex.com
|
|||
//
|
|||
// Copyright (c) 2009-2010 Math.NET
|
|||
//
|
|||
// Permission is hereby granted, free of charge, to any person
|
|||
// obtaining a copy of this software and associated documentation
|
|||
// files (the "Software"), to deal in the Software without
|
|||
// restriction, including without limitation the rights to use,
|
|||
// copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|||
// copies of the Software, and to permit persons to whom the
|
|||
// Software is furnished to do so, subject to the following
|
|||
// conditions:
|
|||
//
|
|||
// The above copyright notice and this permission notice shall be
|
|||
// included in all copies or substantial portions of the Software.
|
|||
//
|
|||
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
|||
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
|
|||
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
|||
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
|
|||
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
|
|||
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
|||
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
|
|||
// OTHER DEALINGS IN THE SOFTWARE.
|
|||
// </copyright>
|
|||
|
|||
namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|||
{ |
|||
using System; |
|||
using System.Linq; |
|||
using Generic; |
|||
using Generic.Factorization; |
|||
using Numerics; |
|||
using Properties; |
|||
|
|||
/// <summary>
|
|||
/// <para>A class which encapsulates the functionality of the singular value decomposition (SVD).</para>
|
|||
/// <para>Suppose M is an m-by-n matrix whose entries are real numbers.
|
|||
/// Then there exists a factorization of the form M = UΣVT where:
|
|||
/// - U is an m-by-m unitary matrix;
|
|||
/// - Σ is m-by-n diagonal matrix with nonnegative real numbers on the diagonal;
|
|||
/// - VT denotes transpose of V, an n-by-n unitary matrix;
|
|||
/// Such a factorization is called a singular-value decomposition of M. A common convention is to order the diagonal
|
|||
/// entries Σ(i,i) in descending order. In this case, the diagonal matrix Σ is uniquely determined
|
|||
/// by M (though the matrices U and V are not). The diagonal entries of Σ are known as the singular values of M.</para>
|
|||
/// </summary>
|
|||
/// <remarks>
|
|||
/// The computation of the singular value decomposition is done at construction time.
|
|||
/// </remarks>
|
|||
public abstract class Svd : Svd<Complex32> |
|||
{ |
|||
/// <summary>
|
|||
/// Gets the effective numerical matrix rank.
|
|||
/// </summary>
|
|||
/// <value>The number of non-negligible singular values.</value>
|
|||
public override int Rank |
|||
{ |
|||
get |
|||
{ |
|||
return VectorS.Count(t => !t.Magnitude.AlmostEqual(0.0f)); |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Gets the two norm of the <see cref="Matrix{T}"/>.
|
|||
/// </summary>
|
|||
/// <returns>The 2-norm of the <see cref="Matrix{T}"/>.</returns>
|
|||
public override Complex32 Norm2 |
|||
{ |
|||
get |
|||
{ |
|||
return VectorS[0].Magnitude; |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Gets the condition number <b>max(S) / min(S)</b>
|
|||
/// </summary>
|
|||
/// <returns>The condition number.</returns>
|
|||
public override Complex32 ConditionNumber |
|||
{ |
|||
get |
|||
{ |
|||
var tmp = Math.Min(MatrixU.RowCount, MatrixVT.ColumnCount) - 1; |
|||
return VectorS[0].Magnitude / VectorS[tmp].Magnitude; |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Gets the determinant of the square matrix for which the SVD was computed.
|
|||
/// </summary>
|
|||
public override Complex32 Determinant |
|||
{ |
|||
get |
|||
{ |
|||
if (MatrixU.RowCount != MatrixVT.ColumnCount) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentMatrixSquare); |
|||
} |
|||
|
|||
var det = Complex32.One; |
|||
foreach (var value in VectorS) |
|||
{ |
|||
det *= value; |
|||
if (value.Magnitude.AlmostEqual(0.0f)) |
|||
{ |
|||
return 0; |
|||
} |
|||
} |
|||
|
|||
return det.Magnitude; |
|||
} |
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,81 @@ |
|||
// <copyright file="Cholesky.cs" company="Math.NET">
|
|||
// Math.NET Numerics, part of the Math.NET Project
|
|||
// http://numerics.mathdotnet.com
|
|||
// http://github.com/mathnet/mathnet-numerics
|
|||
// http://mathnetnumerics.codeplex.com
|
|||
//
|
|||
// Copyright (c) 2009-2010 Math.NET
|
|||
//
|
|||
// Permission is hereby granted, free of charge, to any person
|
|||
// obtaining a copy of this software and associated documentation
|
|||
// files (the "Software"), to deal in the Software without
|
|||
// restriction, including without limitation the rights to use,
|
|||
// copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|||
// copies of the Software, and to permit persons to whom the
|
|||
// Software is furnished to do so, subject to the following
|
|||
// conditions:
|
|||
//
|
|||
// The above copyright notice and this permission notice shall be
|
|||
// included in all copies or substantial portions of the Software.
|
|||
//
|
|||
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
|||
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
|
|||
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
|||
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
|
|||
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
|
|||
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
|||
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
|
|||
// OTHER DEALINGS IN THE SOFTWARE.
|
|||
// </copyright>
|
|||
|
|||
namespace MathNet.Numerics.LinearAlgebra.Double.Factorization |
|||
{ |
|||
using System; |
|||
using Generic.Factorization; |
|||
|
|||
/// <summary>
|
|||
/// <para>A class which encapsulates the functionality of a Cholesky factorization.</para>
|
|||
/// <para>For a symmetric, positive definite matrix A, the Cholesky factorization
|
|||
/// is an lower triangular matrix L so that A = L*L'.</para>
|
|||
/// </summary>
|
|||
/// <remarks>
|
|||
/// The computation of the Cholesky factorization is done at construction time. If the matrix is not symmetric
|
|||
/// or positive definite, the constructor will throw an exception.
|
|||
/// </remarks>
|
|||
public abstract class Cholesky : Cholesky<double> |
|||
{ |
|||
/// <summary>
|
|||
/// Gets the determinant of the matrix for which the Cholesky matrix was computed.
|
|||
/// </summary>
|
|||
public override double Determinant |
|||
{ |
|||
get |
|||
{ |
|||
var det = 1.0; |
|||
for (var j = 0; j < CholeskyFactor.RowCount; j++) |
|||
{ |
|||
det *= CholeskyFactor[j, j] * CholeskyFactor[j, j]; |
|||
} |
|||
|
|||
return det; |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Gets the log determinant of the matrix for which the Cholesky matrix was computed.
|
|||
/// </summary>
|
|||
public override double DeterminantLn |
|||
{ |
|||
get |
|||
{ |
|||
var det = 0.0; |
|||
for (var j = 0; j < CholeskyFactor.RowCount; j++) |
|||
{ |
|||
det += 2 * Math.Log(CholeskyFactor[j, j]); |
|||
} |
|||
|
|||
return det; |
|||
} |
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,114 @@ |
|||
// <copyright file="Evd.cs" company="Math.NET">
|
|||
// Math.NET Numerics, part of the Math.NET Project
|
|||
// http://numerics.mathdotnet.com
|
|||
// http://github.com/mathnet/mathnet-numerics
|
|||
// http://mathnetnumerics.codeplex.com
|
|||
// Copyright (c) 2009-2010 Math.NET
|
|||
// Permission is hereby granted, free of charge, to any person
|
|||
// obtaining a copy of this software and associated documentation
|
|||
// files (the "Software"), to deal in the Software without
|
|||
// restriction, including without limitation the rights to use,
|
|||
// copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|||
// copies of the Software, and to permit persons to whom the
|
|||
// Software is furnished to do so, subject to the following
|
|||
// conditions:
|
|||
// The above copyright notice and this permission notice shall be
|
|||
// included in all copies or substantial portions of the Software.
|
|||
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
|||
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
|
|||
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
|||
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
|
|||
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
|
|||
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
|||
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
|
|||
// OTHER DEALINGS IN THE SOFTWARE.
|
|||
// </copyright>
|
|||
|
|||
namespace MathNet.Numerics.LinearAlgebra.Double.Factorization |
|||
{ |
|||
using System.Numerics; |
|||
using Generic.Factorization; |
|||
|
|||
/// <summary>
|
|||
/// Eigenvalues and eigenvectors of a real matrix.
|
|||
/// </summary>
|
|||
/// <remarks>
|
|||
/// If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is
|
|||
/// diagonal and the eigenvector matrix V is orthogonal.
|
|||
/// I.e. A = V*D*V' and V*VT=I.
|
|||
/// If A is not symmetric, then the eigenvalue matrix D is block diagonal
|
|||
/// with the real eigenvalues in 1-by-1 blocks and any complex eigenvalues,
|
|||
/// lambda + i*mu, in 2-by-2 blocks, [lambda, mu; -mu, lambda]. The
|
|||
/// columns of V represent the eigenvectors in the sense that A*V = V*D,
|
|||
/// i.e. A.Multiply(V) equals V.Multiply(D). The matrix V may be badly
|
|||
/// conditioned, or even singular, so the validity of the equation
|
|||
/// A = V*D*Inverse(V) depends upon V.Condition().
|
|||
/// </remarks>
|
|||
public abstract class Evd : Evd<double> |
|||
{ |
|||
/// <summary>
|
|||
/// Gets the absolute value of determinant of the square matrix for which the EVD was computed.
|
|||
/// </summary>
|
|||
public override double Determinant |
|||
{ |
|||
get |
|||
{ |
|||
var det = Complex.One; |
|||
for (var i = 0; i < VectorEv.Count; i++) |
|||
{ |
|||
det *= VectorEv[i]; |
|||
|
|||
if (VectorEv[i].AlmostEqual(Complex.Zero)) |
|||
{ |
|||
return 0; |
|||
} |
|||
} |
|||
|
|||
return det.Magnitude; |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Gets the effective numerical matrix rank.
|
|||
/// </summary>
|
|||
/// <value>The number of non-negligible singular values.</value>
|
|||
public override int Rank |
|||
{ |
|||
get |
|||
{ |
|||
var rank = 0; |
|||
for (var i = 0; i < VectorEv.Count; i++) |
|||
{ |
|||
if (VectorEv[i].AlmostEqual(Complex.Zero)) |
|||
{ |
|||
continue; |
|||
} |
|||
|
|||
rank++; |
|||
} |
|||
|
|||
return rank; |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Gets a value indicating whether the matrix is full rank or not.
|
|||
/// </summary>
|
|||
/// <value><c>true</c> if the matrix is full rank; otherwise <c>false</c>.</value>
|
|||
public override bool IsFullRank |
|||
{ |
|||
get |
|||
{ |
|||
for (var i = 0; i < VectorEv.Count; i++) |
|||
{ |
|||
if (VectorEv[i].AlmostEqual(Complex.Zero)) |
|||
{ |
|||
return false; |
|||
} |
|||
} |
|||
|
|||
return true; |
|||
} |
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,88 @@ |
|||
// <copyright file="GramSchmidt.cs" company="Math.NET">
|
|||
// Math.NET Numerics, part of the Math.NET Project
|
|||
// http://numerics.mathdotnet.com
|
|||
// http://github.com/mathnet/mathnet-numerics
|
|||
// http://mathnetnumerics.codeplex.com
|
|||
// Copyright (c) 2009-2010 Math.NET
|
|||
// Permission is hereby granted, free of charge, to any person
|
|||
// obtaining a copy of this software and associated documentation
|
|||
// files (the "Software"), to deal in the Software without
|
|||
// restriction, including without limitation the rights to use,
|
|||
// copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|||
// copies of the Software, and to permit persons to whom the
|
|||
// Software is furnished to do so, subject to the following
|
|||
// conditions:
|
|||
// The above copyright notice and this permission notice shall be
|
|||
// included in all copies or substantial portions of the Software.
|
|||
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
|||
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
|
|||
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
|||
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
|
|||
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
|
|||
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
|||
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
|
|||
// OTHER DEALINGS IN THE SOFTWARE.
|
|||
// </copyright>
|
|||
|
|||
namespace MathNet.Numerics.LinearAlgebra.Double.Factorization |
|||
{ |
|||
using System; |
|||
using Generic.Factorization; |
|||
using Properties; |
|||
|
|||
/// <summary>
|
|||
/// <para>A class which encapsulates the functionality of the QR decomposition Modified Gram-Schmidt Orthogonalization.</para>
|
|||
/// <para>Any real square matrix A may be decomposed as A = QR where Q is an orthogonal mxn matrix and R is an nxn upper triangular matrix.</para>
|
|||
/// </summary>
|
|||
/// <remarks>
|
|||
/// The computation of the QR decomposition is done at construction time by modified Gram-Schmidt Orthogonalization.
|
|||
/// </remarks>
|
|||
public abstract class GramSchmidt : GramSchmidt<double> |
|||
{ |
|||
/// <summary>
|
|||
/// Gets the absolute determinant value of the matrix for which the QR matrix was computed.
|
|||
/// </summary>
|
|||
public override double Determinant |
|||
{ |
|||
get |
|||
{ |
|||
if (MatrixR.RowCount != MatrixR.ColumnCount) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentMatrixSquare); |
|||
} |
|||
|
|||
var det = 1.0; |
|||
for (var i = 0; i < MatrixR.ColumnCount; i++) |
|||
{ |
|||
det *= MatrixR.At(i, i); |
|||
if (Math.Abs(MatrixR.At(i, i)).AlmostEqual(0.0)) |
|||
{ |
|||
return 0; |
|||
} |
|||
} |
|||
|
|||
return Convert.ToSingle(Math.Abs(det)); |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Gets a value indicating whether the matrix is full rank or not.
|
|||
/// </summary>
|
|||
/// <value><c>true</c> if the matrix is full rank; otherwise <c>false</c>.</value>
|
|||
public override bool IsFullRank |
|||
{ |
|||
get |
|||
{ |
|||
for (var i = 0; i < MatrixR.ColumnCount; i++) |
|||
{ |
|||
if (Math.Abs(MatrixR.At(i, i)).AlmostEqual(0.0)) |
|||
{ |
|||
return false; |
|||
} |
|||
} |
|||
|
|||
return true; |
|||
} |
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,67 @@ |
|||
// <copyright file="LU.cs" company="Math.NET">
|
|||
// Math.NET Numerics, part of the Math.NET Project
|
|||
// http://numerics.mathdotnet.com
|
|||
// http://github.com/mathnet/mathnet-numerics
|
|||
// http://mathnetnumerics.codeplex.com
|
|||
// Copyright (c) 2009-2010 Math.NET
|
|||
// Permission is hereby granted, free of charge, to any person
|
|||
// obtaining a copy of this software and associated documentation
|
|||
// files (the "Software"), to deal in the Software without
|
|||
// restriction, including without limitation the rights to use,
|
|||
// copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|||
// copies of the Software, and to permit persons to whom the
|
|||
// Software is furnished to do so, subject to the following
|
|||
// conditions:
|
|||
// The above copyright notice and this permission notice shall be
|
|||
// included in all copies or substantial portions of the Software.
|
|||
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
|||
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
|
|||
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
|||
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
|
|||
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
|
|||
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
|||
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
|
|||
// OTHER DEALINGS IN THE SOFTWARE.
|
|||
// </copyright>
|
|||
|
|||
namespace MathNet.Numerics.LinearAlgebra.Double.Factorization |
|||
{ |
|||
using Generic.Factorization; |
|||
|
|||
/// <summary>
|
|||
/// <para>A class which encapsulates the functionality of an LU factorization.</para>
|
|||
/// <para>For a matrix A, the LU factorization is a pair of lower triangular matrix L and
|
|||
/// upper triangular matrix U so that A = L*U.</para>
|
|||
/// <para>In the Math.Net implementation we also store a set of pivot elements for increased
|
|||
/// numerical stability. The pivot elements encode a permutation matrix P such that P*A = L*U.</para>
|
|||
/// </summary>
|
|||
/// <remarks>
|
|||
/// The computation of the LU factorization is done at construction time.
|
|||
/// </remarks>
|
|||
public abstract class LU : LU<double> |
|||
{ |
|||
/// <summary>
|
|||
/// Gets the determinant of the matrix for which the LU factorization was computed.
|
|||
/// </summary>
|
|||
public override double Determinant |
|||
{ |
|||
get |
|||
{ |
|||
var det = 1.0; |
|||
for (var j = 0; j < Factors.RowCount; j++) |
|||
{ |
|||
if (Pivots[j] != j) |
|||
{ |
|||
det *= -Factors.At(j, j); |
|||
} |
|||
else |
|||
{ |
|||
det *= Factors.At(j, j); |
|||
} |
|||
} |
|||
|
|||
return det; |
|||
} |
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,90 @@ |
|||
// <copyright file="QR.cs" company="Math.NET">
|
|||
// Math.NET Numerics, part of the Math.NET Project
|
|||
// http://numerics.mathdotnet.com
|
|||
// http://github.com/mathnet/mathnet-numerics
|
|||
// http://mathnetnumerics.codeplex.com
|
|||
// Copyright (c) 2009-2010 Math.NET
|
|||
// Permission is hereby granted, free of charge, to any person
|
|||
// obtaining a copy of this software and associated documentation
|
|||
// files (the "Software"), to deal in the Software without
|
|||
// restriction, including without limitation the rights to use,
|
|||
// copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|||
// copies of the Software, and to permit persons to whom the
|
|||
// Software is furnished to do so, subject to the following
|
|||
// conditions:
|
|||
// The above copyright notice and this permission notice shall be
|
|||
// included in all copies or substantial portions of the Software.
|
|||
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
|||
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
|
|||
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
|||
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
|
|||
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
|
|||
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
|||
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
|
|||
// OTHER DEALINGS IN THE SOFTWARE.
|
|||
// </copyright>
|
|||
|
|||
namespace MathNet.Numerics.LinearAlgebra.Double.Factorization |
|||
{ |
|||
using System; |
|||
using Generic.Factorization; |
|||
using Properties; |
|||
|
|||
/// <summary>
|
|||
/// <para>A class which encapsulates the functionality of the QR decomposition.</para>
|
|||
/// <para>Any real square matrix A (m x n) may be decomposed as A = QR where Q is an orthogonal matrix (m x m)
|
|||
/// (its columns are orthogonal unit vectors meaning QTQ = I) and R (m x n) is an upper triangular matrix
|
|||
/// (also called right triangular matrix).</para>
|
|||
/// </summary>
|
|||
/// <remarks>
|
|||
/// The computation of the QR decomposition is done at construction time by Householder transformation.
|
|||
/// </remarks>
|
|||
public abstract class QR : QR<double> |
|||
{ |
|||
/// <summary>
|
|||
/// Gets the absolute determinant value of the matrix for which the QR matrix was computed.
|
|||
/// </summary>
|
|||
public override double Determinant |
|||
{ |
|||
get |
|||
{ |
|||
if (MatrixR.RowCount != MatrixR.ColumnCount) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentMatrixSquare); |
|||
} |
|||
|
|||
var det = 1.0; |
|||
for (var i = 0; i < MatrixR.ColumnCount; i++) |
|||
{ |
|||
det *= MatrixR.At(i, i); |
|||
if (Math.Abs(MatrixR.At(i, i)).AlmostEqual(0.0)) |
|||
{ |
|||
return 0; |
|||
} |
|||
} |
|||
|
|||
return Math.Abs(det); |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Gets a value indicating whether the matrix is full rank or not.
|
|||
/// </summary>
|
|||
/// <value><c>true</c> if the matrix is full rank; otherwise <c>false</c>.</value>
|
|||
public override bool IsFullRank |
|||
{ |
|||
get |
|||
{ |
|||
for (var i = 0; i < MatrixR.ColumnCount; i++) |
|||
{ |
|||
if (Math.Abs(MatrixR.At(i, i)).AlmostEqual(0.0)) |
|||
{ |
|||
return false; |
|||
} |
|||
} |
|||
|
|||
return true; |
|||
} |
|||
} |
|||
} |
|||
} |
|||
@ -1,258 +0,0 @@ |
|||
// <copyright file="SparseCholesky.cs" company="Math.NET">
|
|||
// Math.NET Numerics, part of the Math.NET Project
|
|||
// http://numerics.mathdotnet.com
|
|||
// http://github.com/mathnet/mathnet-numerics
|
|||
// http://mathnetnumerics.codeplex.com
|
|||
//
|
|||
// Copyright (c) 2009-2010 Math.NET
|
|||
//
|
|||
// Permission is hereby granted, free of charge, to any person
|
|||
// obtaining a copy of this software and associated documentation
|
|||
// files (the "Software"), to deal in the Software without
|
|||
// restriction, including without limitation the rights to use,
|
|||
// copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|||
// copies of the Software, and to permit persons to whom the
|
|||
// Software is furnished to do so, subject to the following
|
|||
// conditions:
|
|||
//
|
|||
// The above copyright notice and this permission notice shall be
|
|||
// included in all copies or substantial portions of the Software.
|
|||
//
|
|||
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
|||
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
|
|||
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
|||
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
|
|||
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
|
|||
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
|||
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
|
|||
// OTHER DEALINGS IN THE SOFTWARE.
|
|||
// </copyright>
|
|||
|
|||
namespace MathNet.Numerics.LinearAlgebra.Double.Factorization |
|||
{ |
|||
using System; |
|||
using Generic; |
|||
using Generic.Factorization; |
|||
using Properties; |
|||
|
|||
/// <summary>
|
|||
/// <para>A class which encapsulates the functionality of a Cholesky factorization for soarse matrices.</para>
|
|||
/// <para>For a symmetric, positive definite matrix A, the Cholesky factorization
|
|||
/// is an lower triangular matrix L so that A = L*L'.</para>
|
|||
/// </summary>
|
|||
/// <remarks>
|
|||
/// The computation of the Cholesky factorization is done at construction time. If the matrix is not symmetric
|
|||
/// or positive definite, the constructor will throw an exception.
|
|||
/// </remarks>
|
|||
public class SparseCholesky : Cholesky<double> |
|||
{ |
|||
/// <summary>
|
|||
/// Initializes a new instance of the <see cref="SparseCholesky"/> class. This object will compute the
|
|||
/// Cholesky factorization when the constructor is called and cache it's factorization.
|
|||
/// </summary>
|
|||
/// <param name="matrix">The matrix to factor.</param>
|
|||
/// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
|
|||
/// <exception cref="ArgumentException">If <paramref name="matrix"/> is not a square matrix.</exception>
|
|||
/// <exception cref="ArgumentException">If <paramref name="matrix"/> is not positive definite.</exception>
|
|||
public SparseCholesky(Matrix<double> matrix) |
|||
{ |
|||
if (matrix == null) |
|||
{ |
|||
throw new ArgumentNullException("matrix"); |
|||
} |
|||
|
|||
if (matrix.RowCount != matrix.ColumnCount) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentMatrixSquare); |
|||
} |
|||
|
|||
// Create a new matrix for the Cholesky factor, then perform factorization (while overwriting).
|
|||
CholeskyFactor = matrix.Clone(); |
|||
for (var j = 0; j < CholeskyFactor.RowCount; j++) |
|||
{ |
|||
var d = 0.0; |
|||
for (var k = 0; k < j; k++) |
|||
{ |
|||
var s = 0.0; |
|||
for (var i = 0; i < k; i++) |
|||
{ |
|||
s += CholeskyFactor.At(k, i) * CholeskyFactor.At(j, i); |
|||
} |
|||
|
|||
s = (matrix.At(j, k) - s) / CholeskyFactor.At(k, k); |
|||
CholeskyFactor.At(j, k, s); |
|||
d += s * s; |
|||
} |
|||
|
|||
d = matrix.At(j, j) - d; |
|||
if (d <= 0.0) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentMatrixPositiveDefinite); |
|||
} |
|||
|
|||
CholeskyFactor.At(j, j, Math.Sqrt(d)); |
|||
for (var k = j + 1; k < CholeskyFactor.RowCount; k++) |
|||
{ |
|||
CholeskyFactor.At(j, k, 0.0); |
|||
} |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Solves a system of linear equations, <b>AX = B</b>, with A Cholesky factorized.
|
|||
/// </summary>
|
|||
/// <param name="input">The right hand side <see cref="Matrix{T}"/>, <b>B</b>.</param>
|
|||
/// <param name="result">The left hand side <see cref="Matrix{T}"/>, <b>X</b>.</param>
|
|||
public override void Solve(Matrix<double> input, Matrix<double> result) |
|||
{ |
|||
if (input == null) |
|||
{ |
|||
throw new ArgumentNullException("input"); |
|||
} |
|||
|
|||
if (result == null) |
|||
{ |
|||
throw new ArgumentNullException("result"); |
|||
} |
|||
|
|||
// Check for proper dimensions.
|
|||
if (result.RowCount != input.RowCount) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension); |
|||
} |
|||
|
|||
if (result.ColumnCount != input.ColumnCount) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension); |
|||
} |
|||
|
|||
if (input.RowCount != CholeskyFactor.RowCount) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentMatrixDimensions); |
|||
} |
|||
|
|||
input.CopyTo(result); |
|||
var order = CholeskyFactor.RowCount; |
|||
|
|||
for (var c = 0; c < result.ColumnCount; c++) |
|||
{ |
|||
// Solve L*Y = B;
|
|||
double sum; |
|||
for (var i = 0; i < order; i++) |
|||
{ |
|||
sum = result.At(i, c); |
|||
for (var k = i - 1; k >= 0; k--) |
|||
{ |
|||
sum -= CholeskyFactor.At(i, k) * result.At(k, c); |
|||
} |
|||
|
|||
result.At(i, c, sum / CholeskyFactor.At(i, i)); |
|||
} |
|||
|
|||
// Solve L'*X = Y;
|
|||
for (var i = order - 1; i >= 0; i--) |
|||
{ |
|||
sum = result.At(i, c); |
|||
for (var k = i + 1; k < order; k++) |
|||
{ |
|||
sum -= CholeskyFactor.At(k, i) * result.At(k, c); |
|||
} |
|||
|
|||
result.At(i, c, sum / CholeskyFactor.At(i, i)); |
|||
} |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Solves a system of linear equations, <b>Ax = b</b>, with A Cholesky factorized.
|
|||
/// </summary>
|
|||
/// <param name="input">The right hand side vector, <b>b</b>.</param>
|
|||
/// <param name="result">The left hand side <see cref="Matrix{T}"/>, <b>x</b>.</param>
|
|||
public override void Solve(Vector<double> input, Vector<double> result) |
|||
{ |
|||
// Check for proper arguments.
|
|||
if (input == null) |
|||
{ |
|||
throw new ArgumentNullException("input"); |
|||
} |
|||
|
|||
if (result == null) |
|||
{ |
|||
throw new ArgumentNullException("result"); |
|||
} |
|||
|
|||
// Check for proper dimensions.
|
|||
if (input.Count != result.Count) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentVectorsSameLength); |
|||
} |
|||
|
|||
if (input.Count != CholeskyFactor.RowCount) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentMatrixDimensions); |
|||
} |
|||
|
|||
input.CopyTo(result); |
|||
var order = CholeskyFactor.RowCount; |
|||
|
|||
// Solve L*Y = B;
|
|||
double sum; |
|||
for (var i = 0; i < order; i++) |
|||
{ |
|||
sum = result[i]; |
|||
for (var k = i - 1; k >= 0; k--) |
|||
{ |
|||
sum -= CholeskyFactor.At(i, k) * result[k]; |
|||
} |
|||
|
|||
result[i] = sum / CholeskyFactor.At(i, i); |
|||
} |
|||
|
|||
// Solve L'*X = Y;
|
|||
for (var i = order - 1; i >= 0; i--) |
|||
{ |
|||
sum = result[i]; |
|||
for (var k = i + 1; k < order; k++) |
|||
{ |
|||
sum -= CholeskyFactor.At(k, i) * result[k]; |
|||
} |
|||
|
|||
result[i] = sum / CholeskyFactor.At(i, i); |
|||
} |
|||
} |
|||
|
|||
#region Simple arithmetic of type T
|
|||
/// <summary>
|
|||
/// Add two values T+T
|
|||
/// </summary>
|
|||
/// <param name="val1">Left operand value</param>
|
|||
/// <param name="val2">Right operand value</param>
|
|||
/// <returns>Result of addition</returns>
|
|||
protected sealed override double AddT(double val1, double val2) |
|||
{ |
|||
return val1 + val2; |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Multiply two values T*T
|
|||
/// </summary>
|
|||
/// <param name="val1">Left operand value</param>
|
|||
/// <param name="val2">Right operand value</param>
|
|||
/// <returns>Result of multiplication</returns>
|
|||
protected sealed override double MultiplyT(double val1, double val2) |
|||
{ |
|||
return val1 * val2; |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Returns the natural (base e) logarithm of a specified number.
|
|||
/// </summary>
|
|||
/// <param name="val1"> A number whose logarithm is to be found</param>
|
|||
/// <returns>Natural (base e) logarithm </returns>
|
|||
protected sealed override double LogT(double val1) |
|||
{ |
|||
return Math.Log(val1); |
|||
} |
|||
#endregion
|
|||
} |
|||
} |
|||
@ -1,317 +0,0 @@ |
|||
// <copyright file="SparseLU.cs" company="Math.NET">
|
|||
// Math.NET Numerics, part of the Math.NET Project
|
|||
// http://numerics.mathdotnet.com
|
|||
// http://github.com/mathnet/mathnet-numerics
|
|||
// http://mathnetnumerics.codeplex.com
|
|||
//
|
|||
// Copyright (c) 2009-2010 Math.NET
|
|||
//
|
|||
// Permission is hereby granted, free of charge, to any person
|
|||
// obtaining a copy of this software and associated documentation
|
|||
// files (the "Software"), to deal in the Software without
|
|||
// restriction, including without limitation the rights to use,
|
|||
// copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|||
// copies of the Software, and to permit persons to whom the
|
|||
// Software is furnished to do so, subject to the following
|
|||
// conditions:
|
|||
//
|
|||
// The above copyright notice and this permission notice shall be
|
|||
// included in all copies or substantial portions of the Software.
|
|||
//
|
|||
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
|||
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
|
|||
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
|||
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
|
|||
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
|
|||
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
|||
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
|
|||
// OTHER DEALINGS IN THE SOFTWARE.
|
|||
// </copyright>
|
|||
|
|||
namespace MathNet.Numerics.LinearAlgebra.Double.Factorization |
|||
{ |
|||
using System; |
|||
using Generic; |
|||
using Generic.Factorization; |
|||
using Properties; |
|||
|
|||
/// <summary>
|
|||
/// <para>A class which encapsulates the functionality of an LU factorization.</para>
|
|||
/// <para>For a matrix A, the LU factorization is a pair of lower triangular matrix L and
|
|||
/// upper triangular matrix U so that A = L*U.</para>
|
|||
/// </summary>
|
|||
/// <remarks>
|
|||
/// The computation of the LU factorization is done at construction time.
|
|||
/// </remarks>
|
|||
public class SparseLU : LU<double> |
|||
{ |
|||
/// <summary>
|
|||
/// Initializes a new instance of the <see cref="SparseLU"/> class. This object will compute the
|
|||
/// LU factorization when the constructor is called and cache it's factorization.
|
|||
/// </summary>
|
|||
/// <param name="matrix">The matrix to factor.</param>
|
|||
/// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
|
|||
/// <exception cref="ArgumentException">If <paramref name="matrix"/> is not a square matrix.</exception>
|
|||
public SparseLU(Matrix<double> matrix) |
|||
{ |
|||
if (matrix == null) |
|||
{ |
|||
throw new ArgumentNullException("matrix"); |
|||
} |
|||
|
|||
if (matrix.RowCount != matrix.ColumnCount) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentMatrixSquare); |
|||
} |
|||
|
|||
// Create an array for the pivot indices.
|
|||
var order = matrix.RowCount; |
|||
Factors = matrix.Clone(); |
|||
Pivots = new int[order]; |
|||
|
|||
// Initialize the pivot matrix to the identity permutation.
|
|||
for (var i = 0; i < order; i++) |
|||
{ |
|||
Pivots[i] = i; |
|||
} |
|||
|
|||
var vectorLUcolj = new double[order]; |
|||
for (var j = 0; j < order; j++) |
|||
{ |
|||
// Make a copy of the j-th column to localize references.
|
|||
for (var i = 0; i < order; i++) |
|||
{ |
|||
vectorLUcolj[i] = Factors.At(i, j); |
|||
} |
|||
|
|||
// Apply previous transformations.
|
|||
for (var i = 0; i < order; i++) |
|||
{ |
|||
var kmax = Math.Min(i, j); |
|||
var s = 0.0; |
|||
for (var k = 0; k < kmax; k++) |
|||
{ |
|||
s += Factors.At(i, k) * vectorLUcolj[k]; |
|||
} |
|||
|
|||
vectorLUcolj[i] -= s; |
|||
Factors.At(i, j, vectorLUcolj[i]); |
|||
} |
|||
|
|||
// Find pivot and exchange if necessary.
|
|||
var p = j; |
|||
for (var i = j + 1; i < order; i++) |
|||
{ |
|||
if (Math.Abs(vectorLUcolj[i]) > Math.Abs(vectorLUcolj[p])) |
|||
{ |
|||
p = i; |
|||
} |
|||
} |
|||
|
|||
if (p != j) |
|||
{ |
|||
for (var k = 0; k < order; k++) |
|||
{ |
|||
var temp = Factors.At(p, k); |
|||
Factors.At(p, k, Factors.At(j, k)); |
|||
Factors.At(j, k, temp); |
|||
} |
|||
|
|||
Pivots[j] = p; |
|||
} |
|||
|
|||
// Compute multipliers.
|
|||
if (j < order & Factors.At(j, j) != 0.0) |
|||
{ |
|||
for (var i = j + 1; i < order; i++) |
|||
{ |
|||
Factors.At(i, j, (Factors.At(i, j) / Factors.At(j, j))); |
|||
} |
|||
} |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Solves a system of linear equations, <c>AX = B</c>, with A LU factorized.
|
|||
/// </summary>
|
|||
/// <param name="input">The right hand side <see cref="Matrix{T}"/>, <c>B</c>.</param>
|
|||
/// <param name="result">The left hand side <see cref="Matrix{T}"/>, <c>X</c>.</param>
|
|||
public override void Solve(Matrix<double> input, Matrix<double> result) |
|||
{ |
|||
// Check for proper arguments.
|
|||
if (input == null) |
|||
{ |
|||
throw new ArgumentNullException("input"); |
|||
} |
|||
|
|||
if (result == null) |
|||
{ |
|||
throw new ArgumentNullException("result"); |
|||
} |
|||
|
|||
// Check for proper dimensions.
|
|||
if (result.RowCount != input.RowCount) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension); |
|||
} |
|||
|
|||
if (result.ColumnCount != input.ColumnCount) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension); |
|||
} |
|||
|
|||
if (input.RowCount != Factors.RowCount) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentMatrixDimensions); |
|||
} |
|||
|
|||
// Copy the contents of input to result.
|
|||
input.CopyTo(result); |
|||
for (var i = 0; i < Pivots.Length; i++) |
|||
{ |
|||
if (Pivots[i] == i) |
|||
{ |
|||
continue; |
|||
} |
|||
|
|||
var p = Pivots[i]; |
|||
for (var j = 0; j < result.ColumnCount; j++) |
|||
{ |
|||
var temp = result.At(p, j); |
|||
result.At(p, j, result.At(i, j)); |
|||
result.At(i, j, temp); |
|||
} |
|||
} |
|||
|
|||
var order = Factors.RowCount; |
|||
|
|||
// Solve L*Y = P*B
|
|||
for (var k = 0; k < order; k++) |
|||
{ |
|||
for (var i = k + 1; i < order; i++) |
|||
{ |
|||
for (var j = 0; j < result.ColumnCount; j++) |
|||
{ |
|||
var temp = result.At(k, j) * Factors.At(i, k); |
|||
result.At(i, j, result.At(i, j) - temp); |
|||
} |
|||
} |
|||
} |
|||
|
|||
// Solve U*X = Y;
|
|||
for (var k = order - 1; k >= 0; k--) |
|||
{ |
|||
for (var j = 0; j < result.ColumnCount; j++) |
|||
{ |
|||
result.At(k, j, (result.At(k, j) / Factors.At(k, k))); |
|||
} |
|||
|
|||
for (var i = 0; i < k; i++) |
|||
{ |
|||
for (var j = 0; j < result.ColumnCount; j++) |
|||
{ |
|||
var temp = result.At(k, j) * Factors.At(i, k); |
|||
result.At(i, j, result.At(i, j) - temp); |
|||
} |
|||
} |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Solves a system of linear equations, <c>Ax = b</c>, with A LU factorized.
|
|||
/// </summary>
|
|||
/// <param name="input">The right hand side vector, <c>b</c>.</param>
|
|||
/// <param name="result">The left hand side <see cref="Matrix{T}"/>, <c>x</c>.</param>
|
|||
public override void Solve(Vector<double> input, Vector<double> result) |
|||
{ |
|||
// Check for proper arguments.
|
|||
if (input == null) |
|||
{ |
|||
throw new ArgumentNullException("input"); |
|||
} |
|||
|
|||
if (result == null) |
|||
{ |
|||
throw new ArgumentNullException("result"); |
|||
} |
|||
|
|||
// Check for proper dimensions.
|
|||
if (input.Count != result.Count) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentVectorsSameLength); |
|||
} |
|||
|
|||
if (input.Count != Factors.RowCount) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentMatrixDimensions); |
|||
} |
|||
|
|||
// Copy the contents of input to result.
|
|||
input.CopyTo(result); |
|||
for (var i = 0; i < Pivots.Length; i++) |
|||
{ |
|||
if (Pivots[i] == i) |
|||
{ |
|||
continue; |
|||
} |
|||
|
|||
var p = Pivots[i]; |
|||
var temp = result[p]; |
|||
result[p] = result[i]; |
|||
result[i] = temp; |
|||
} |
|||
|
|||
var order = Factors.RowCount; |
|||
|
|||
// Solve L*Y = P*B
|
|||
for (var k = 0; k < order; k++) |
|||
{ |
|||
for (var i = k + 1; i < order; i++) |
|||
{ |
|||
result[i] -= result[k] * Factors.At(i, k); |
|||
} |
|||
} |
|||
|
|||
// Solve U*X = Y;
|
|||
for (var k = order - 1; k >= 0; k--) |
|||
{ |
|||
result[k] /= Factors.At(k, k); |
|||
for (var i = 0; i < k; i++) |
|||
{ |
|||
result[i] -= result[k] * Factors.At(i, k); |
|||
} |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Returns the inverse of this matrix. The inverse is calculated using LU decomposition.
|
|||
/// </summary>
|
|||
/// <returns>The inverse of this matrix.</returns>
|
|||
public override Matrix<double> Inverse() |
|||
{ |
|||
var order = Factors.RowCount; |
|||
var inverse = Factors.CreateMatrix(order, order); |
|||
for (var i = 0; i < order; i++) |
|||
{ |
|||
inverse.At(i, i, 1.0); |
|||
} |
|||
|
|||
return Solve(inverse); |
|||
} |
|||
|
|||
#region Simple arithmetic of type T
|
|||
|
|||
/// <summary>
|
|||
/// Multiply two values T*T
|
|||
/// </summary>
|
|||
/// <param name="val1">Left operand value</param>
|
|||
/// <param name="val2">Right operand value</param>
|
|||
/// <returns>Result of multiplication</returns>
|
|||
protected sealed override double MultiplyT(double val1, double val2) |
|||
{ |
|||
return val1 * val2; |
|||
} |
|||
|
|||
#endregion
|
|||
} |
|||
} |
|||
@ -1,358 +0,0 @@ |
|||
// <copyright file="SparseQR.cs" company="Math.NET">
|
|||
// Math.NET Numerics, part of the Math.NET Project
|
|||
// http://numerics.mathdotnet.com
|
|||
// http://github.com/mathnet/mathnet-numerics
|
|||
// http://mathnetnumerics.codeplex.com
|
|||
//
|
|||
// Copyright (c) 2009-2010 Math.NET
|
|||
//
|
|||
// Permission is hereby granted, free of charge, to any person
|
|||
// obtaining a copy of this software and associated documentation
|
|||
// files (the "Software"), to deal in the Software without
|
|||
// restriction, including without limitation the rights to use,
|
|||
// copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|||
// copies of the Software, and to permit persons to whom the
|
|||
// Software is furnished to do so, subject to the following
|
|||
// conditions:
|
|||
//
|
|||
// The above copyright notice and this permission notice shall be
|
|||
// included in all copies or substantial portions of the Software.
|
|||
//
|
|||
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
|||
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
|
|||
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
|||
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
|
|||
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
|
|||
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
|||
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
|
|||
// OTHER DEALINGS IN THE SOFTWARE.
|
|||
// </copyright>
|
|||
|
|||
namespace MathNet.Numerics.LinearAlgebra.Double.Factorization |
|||
{ |
|||
using System; |
|||
using System.Linq; |
|||
using Generic; |
|||
using Generic.Factorization; |
|||
using Properties; |
|||
|
|||
/// <summary>
|
|||
/// <para>A class which encapsulates the functionality of the QR decomposition.</para>
|
|||
/// <para>Any real square matrix A may be decomposed as A = QR where Q is an orthogonal matrix
|
|||
/// (its columns are orthogonal unit vectors meaning QTQ = I) and R is an upper triangular matrix
|
|||
/// (also called right triangular matrix).</para>
|
|||
/// </summary>
|
|||
/// <remarks>
|
|||
/// The computation of the QR decomposition is done at construction time by Householder transformation.
|
|||
/// </remarks>
|
|||
public class SparseQR : QR<double> |
|||
{ |
|||
/// <summary>
|
|||
/// Initializes a new instance of the <see cref="SparseQR"/> class. This object will compute the
|
|||
/// QR factorization when the constructor is called and cache it's factorization.
|
|||
/// </summary>
|
|||
/// <param name="matrix">The matrix to factor.</param>
|
|||
/// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
|
|||
public SparseQR(Matrix<double> matrix) |
|||
{ |
|||
if (matrix == null) |
|||
{ |
|||
throw new ArgumentNullException("matrix"); |
|||
} |
|||
|
|||
if (matrix.RowCount < matrix.ColumnCount) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentMatrixDimensions); |
|||
} |
|||
|
|||
MatrixR = matrix.Clone(); |
|||
MatrixQ = matrix.CreateMatrix(matrix.RowCount, matrix.RowCount); |
|||
|
|||
for (var i = 0; i < matrix.RowCount; i++) |
|||
{ |
|||
MatrixQ.At(i, i, 1.0); |
|||
} |
|||
|
|||
var minmn = Math.Min(matrix.RowCount, matrix.ColumnCount); |
|||
var u = new double[minmn][]; |
|||
for (var i = 0; i < minmn; i++) |
|||
{ |
|||
u[i] = GenerateColumn(MatrixR, i, matrix.RowCount - 1, i); |
|||
ComputeQR(u[i], MatrixR, i, matrix.RowCount - 1, i + 1, matrix.ColumnCount - 1); |
|||
} |
|||
|
|||
for (var i = minmn - 1; i >= 0; i--) |
|||
{ |
|||
ComputeQR(u[i], MatrixQ, i, matrix.RowCount - 1, i, matrix.RowCount - 1); |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Generate column from initial matrix to work array
|
|||
/// </summary>
|
|||
/// <param name="a">Initial matrix</param>
|
|||
/// <param name="rowStart">The firts row</param>
|
|||
/// <param name="rowEnd">The last row</param>
|
|||
/// <param name="column">Column index</param>
|
|||
/// <returns>Generated vector</returns>
|
|||
private static double[] GenerateColumn(Matrix<double> a, int rowStart, int rowEnd, int column) |
|||
{ |
|||
var ru = rowEnd - rowStart + 1; |
|||
var u = new double[ru]; |
|||
|
|||
for (var i = rowStart; i <= rowEnd; i++) |
|||
{ |
|||
u[i - rowStart] = a.At(i, rowStart); |
|||
a.At(i, rowStart, 0.0); |
|||
} |
|||
|
|||
var norm = u.Sum(t => t * t); |
|||
norm = Math.Sqrt(norm); |
|||
|
|||
if (rowStart == rowEnd || norm == 0) |
|||
{ |
|||
a.At(rowStart, column, -u[0]); |
|||
u[0] = Math.Sqrt(2.0); |
|||
return u; |
|||
} |
|||
|
|||
var scale = 1.0 / norm; |
|||
if (u[0] < 0.0) |
|||
{ |
|||
scale *= -1.0; |
|||
} |
|||
|
|||
a.At(rowStart, column, -1.0 / scale); |
|||
|
|||
for (var i = 0; i < ru; i++) |
|||
{ |
|||
u[i] *= scale; |
|||
} |
|||
|
|||
u[0] += 1.0; |
|||
var s = Math.Sqrt(1.0 / u[0]); |
|||
|
|||
for (var i = 0; i < ru; i++) |
|||
{ |
|||
u[i] *= s; |
|||
} |
|||
|
|||
return u; |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Perform calculation of Q or R
|
|||
/// </summary>
|
|||
/// <param name="u">Work array</param>
|
|||
/// <param name="a">Q or R matrices</param>
|
|||
/// <param name="rowStart">The first row</param>
|
|||
/// <param name="rowEnd">The last row</param>
|
|||
/// <param name="columnStart">The first column</param>
|
|||
/// <param name="columnEnd">The last column</param>
|
|||
private static void ComputeQR(double[] u, Matrix<double> a, int rowStart, int rowEnd, int columnStart, int columnEnd) |
|||
{ |
|||
if (rowEnd < rowStart || columnEnd < columnStart) |
|||
{ |
|||
return; |
|||
} |
|||
|
|||
var v = new double[columnEnd - columnStart + 1]; |
|||
for (var j = columnStart; j <= columnEnd; j++) |
|||
{ |
|||
v[j - columnStart] = 0.0; |
|||
} |
|||
|
|||
for (var i = rowStart; i <= rowEnd; i++) |
|||
{ |
|||
for (var j = columnStart; j <= columnEnd; j++) |
|||
{ |
|||
v[j - columnStart] = v[j - columnStart] + (u[i - rowStart] * a.At(i, j)); |
|||
} |
|||
} |
|||
|
|||
for (var i = rowStart; i <= rowEnd; i++) |
|||
{ |
|||
for (var j = columnStart; j <= columnEnd; j++) |
|||
{ |
|||
a.At(i, j, a.At(i, j) - (u[i - rowStart] * v[j - columnStart])); |
|||
} |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Solves a system of linear equations, <b>AX = B</b>, with A QR factorized.
|
|||
/// </summary>
|
|||
/// <param name="input">The right hand side <see cref="Matrix{T}"/>, <b>B</b>.</param>
|
|||
/// <param name="result">The left hand side <see cref="Matrix{T}"/>, <b>X</b>.</param>
|
|||
public override void Solve(Matrix<double> input, Matrix<double> result) |
|||
{ |
|||
// Check for proper arguments.
|
|||
if (input == null) |
|||
{ |
|||
throw new ArgumentNullException("input"); |
|||
} |
|||
|
|||
if (result == null) |
|||
{ |
|||
throw new ArgumentNullException("result"); |
|||
} |
|||
|
|||
// The solution X should have the same number of columns as B
|
|||
if (input.ColumnCount != result.ColumnCount) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension); |
|||
} |
|||
|
|||
// The dimension compatibility conditions for X = A\B require the two matrices A and B to have the same number of rows
|
|||
if (MatrixR.RowCount != input.RowCount) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension); |
|||
} |
|||
|
|||
// The solution X row dimension is equal to the column dimension of A
|
|||
if (MatrixR.ColumnCount != result.RowCount) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension); |
|||
} |
|||
|
|||
var inputCopy = input.Clone(); |
|||
|
|||
// Compute Y = transpose(Q)*B
|
|||
var bn = inputCopy.ColumnCount; |
|||
var column = new double[MatrixR.RowCount]; |
|||
for (var j = 0; j < bn; j++) |
|||
{ |
|||
for (var k = 0; k < MatrixR.RowCount; k++) |
|||
{ |
|||
column[k] = inputCopy.At(k, j); |
|||
} |
|||
|
|||
for (var i = 0; i < MatrixR.RowCount; i++) |
|||
{ |
|||
double s = 0; |
|||
for (var k = 0; k < MatrixR.RowCount; k++) |
|||
{ |
|||
s += MatrixQ.At(k, i) * column[k]; |
|||
} |
|||
|
|||
inputCopy.At(i, j, s); |
|||
} |
|||
} |
|||
|
|||
// Solve R*X = Y;
|
|||
for (var k = MatrixR.ColumnCount - 1; k >= 0; k--) |
|||
{ |
|||
for (var j = 0; j < bn; j++) |
|||
{ |
|||
inputCopy.At(k, j, inputCopy.At(k, j) / MatrixR.At(k, k)); |
|||
} |
|||
|
|||
for (var i = 0; i < k; i++) |
|||
{ |
|||
for (var j = 0; j < bn; j++) |
|||
{ |
|||
inputCopy.At(i, j, inputCopy.At(i, j) - (inputCopy.At(k, j) * MatrixR.At(i, k))); |
|||
} |
|||
} |
|||
} |
|||
|
|||
for (var i = 0; i < MatrixR.ColumnCount; i++) |
|||
{ |
|||
for (var j = 0; j < inputCopy.ColumnCount; j++) |
|||
{ |
|||
result.At(i, j, inputCopy.At(i, j)); |
|||
} |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Solves a system of linear equations, <b>Ax = b</b>, with A QR factorized.
|
|||
/// </summary>
|
|||
/// <param name="input">The right hand side vector, <b>b</b>.</param>
|
|||
/// <param name="result">The left hand side <see cref="Matrix{T}"/>, <b>x</b>.</param>
|
|||
public override void Solve(Vector<double> input, Vector<double> result) |
|||
{ |
|||
if (input == null) |
|||
{ |
|||
throw new ArgumentNullException("input"); |
|||
} |
|||
|
|||
if (result == null) |
|||
{ |
|||
throw new ArgumentNullException("result"); |
|||
} |
|||
|
|||
// Ax=b where A is an m x n matrix
|
|||
// Check that b is a column vector with m entries
|
|||
if (MatrixR.RowCount != input.Count) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentVectorsSameLength); |
|||
} |
|||
|
|||
// Check that x is a column vector with n entries
|
|||
if (MatrixR.ColumnCount != result.Count) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentMatrixDimensions); |
|||
} |
|||
|
|||
var inputCopy = input.Clone(); |
|||
|
|||
// Compute Y = transpose(Q)*B
|
|||
var column = new double[MatrixR.RowCount]; |
|||
for (var k = 0; k < MatrixR.RowCount; k++) |
|||
{ |
|||
column[k] = inputCopy[k]; |
|||
} |
|||
|
|||
for (var i = 0; i < MatrixR.RowCount; i++) |
|||
{ |
|||
double s = 0; |
|||
for (var k = 0; k < MatrixR.RowCount; k++) |
|||
{ |
|||
s += MatrixQ.At(k, i) * column[k]; |
|||
} |
|||
|
|||
inputCopy[i] = s; |
|||
} |
|||
|
|||
// Solve R*X = Y;
|
|||
for (var k = MatrixR.ColumnCount - 1; k >= 0; k--) |
|||
{ |
|||
inputCopy[k] /= MatrixR.At(k, k); |
|||
for (var i = 0; i < k; i++) |
|||
{ |
|||
inputCopy[i] -= inputCopy[k] * MatrixR.At(i, k); |
|||
} |
|||
} |
|||
|
|||
for (var i = 0; i < MatrixR.ColumnCount; i++) |
|||
{ |
|||
result[i] = inputCopy[i]; |
|||
} |
|||
} |
|||
|
|||
#region Simple arithmetic of type T
|
|||
|
|||
/// <summary>
|
|||
/// Multiply two values T*T
|
|||
/// </summary>
|
|||
/// <param name="val1">Left operand value</param>
|
|||
/// <param name="val2">Right operand value</param>
|
|||
/// <returns>Result of multiplication</returns>
|
|||
protected sealed override double MultiplyT(double val1, double val2) |
|||
{ |
|||
return val1 * val2; |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Returns the absolute value of a specified number.
|
|||
/// </summary>
|
|||
/// <param name="val1"> A number whose absolute is to be found</param>
|
|||
/// <returns>Absolute value </returns>
|
|||
protected sealed override double AbsoluteT(double val1) |
|||
{ |
|||
return Math.Abs(val1); |
|||
} |
|||
#endregion
|
|||
} |
|||
} |
|||
@ -1,950 +0,0 @@ |
|||
// <copyright file="SparseSvd.cs" company="Math.NET">
|
|||
// Math.NET Numerics, part of the Math.NET Project
|
|||
// http://numerics.mathdotnet.com
|
|||
// http://github.com/mathnet/mathnet-numerics
|
|||
// http://mathnetnumerics.codeplex.com
|
|||
//
|
|||
// Copyright (c) 2009-2010 Math.NET
|
|||
//
|
|||
// Permission is hereby granted, free of charge, to any person
|
|||
// obtaining a copy of this software and associated documentation
|
|||
// files (the "Software"), to deal in the Software without
|
|||
// restriction, including without limitation the rights to use,
|
|||
// copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|||
// copies of the Software, and to permit persons to whom the
|
|||
// Software is furnished to do so, subject to the following
|
|||
// conditions:
|
|||
//
|
|||
// The above copyright notice and this permission notice shall be
|
|||
// included in all copies or substantial portions of the Software.
|
|||
//
|
|||
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
|||
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
|
|||
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
|||
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
|
|||
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
|
|||
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
|||
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
|
|||
// OTHER DEALINGS IN THE SOFTWARE.
|
|||
// </copyright>
|
|||
namespace MathNet.Numerics.LinearAlgebra.Double.Factorization |
|||
{ |
|||
using System; |
|||
using Generic; |
|||
using Generic.Factorization; |
|||
using Properties; |
|||
|
|||
/// <summary>
|
|||
/// <para>A class which encapsulates the functionality of the singular value decomposition (SVD) for <see cref="Matrix{T}"/>.</para>
|
|||
/// <para>Suppose M is an m-by-n matrix whose entries are real numbers.
|
|||
/// Then there exists a factorization of the form M = UΣVT where:
|
|||
/// - U is an m-by-m unitary matrix;
|
|||
/// - Σ is m-by-n diagonal matrix with nonnegative real numbers on the diagonal;
|
|||
/// - VT denotes transpose of V, an n-by-n unitary matrix;
|
|||
/// Such a factorization is called a singular-value decomposition of M. A common convention is to order the diagonal
|
|||
/// entries Σ(i,i) in descending order. In this case, the diagonal matrix Σ is uniquely determined
|
|||
/// by M (though the matrices U and V are not). The diagonal entries of Σ are known as the singular values of M.</para>
|
|||
/// </summary>
|
|||
/// <remarks>
|
|||
/// The computation of the singular value decomposition is done at construction time.
|
|||
/// </remarks>
|
|||
public class SparseSvd : Svd<double> |
|||
{ |
|||
/// <summary>
|
|||
/// Initializes a new instance of the <see cref="SparseSvd"/> class. This object will compute the
|
|||
/// the singular value decomposition when the constructor is called and cache it's decomposition.
|
|||
/// </summary>
|
|||
/// <param name="matrix">The matrix to factor.</param>
|
|||
/// <param name="computeVectors">Compute the singular U and VT vectors or not.</param>
|
|||
/// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <b>null</b>.</exception>
|
|||
/// <exception cref="ArgumentException">If SVD algorithm failed to converge with matrix <paramref name="matrix"/>.</exception>
|
|||
public SparseSvd(Matrix<double> matrix, bool computeVectors) |
|||
{ |
|||
if (matrix == null) |
|||
{ |
|||
throw new ArgumentNullException("matrix"); |
|||
} |
|||
|
|||
ComputeVectors = computeVectors; |
|||
var nm = Math.Min(matrix.RowCount + 1, matrix.ColumnCount); |
|||
var matrixCopy = matrix.Clone(); |
|||
|
|||
VectorS = matrixCopy.CreateVector(nm); |
|||
MatrixU = matrixCopy.CreateMatrix(matrixCopy.RowCount, matrixCopy.RowCount); |
|||
MatrixVT = matrixCopy.CreateMatrix(matrixCopy.ColumnCount, matrixCopy.ColumnCount); |
|||
|
|||
const int Maxiter = 1000; |
|||
var e = new double[matrixCopy.ColumnCount]; |
|||
var work = new double[matrixCopy.RowCount]; |
|||
|
|||
int i, j; |
|||
int l, lp1; |
|||
var cs = 0.0; |
|||
var sn = 0.0; |
|||
double t; |
|||
|
|||
var ncu = matrixCopy.RowCount; |
|||
|
|||
// Reduce matrixCopy to bidiagonal form, storing the diagonal elements
|
|||
// In s and the super-diagonal elements in e.
|
|||
var nct = Math.Min(matrixCopy.RowCount - 1, matrixCopy.ColumnCount); |
|||
var nrt = Math.Max(0, Math.Min(matrixCopy.ColumnCount - 2, matrixCopy.RowCount)); |
|||
var lu = Math.Max(nct, nrt); |
|||
for (l = 0; l < lu; l++) |
|||
{ |
|||
lp1 = l + 1; |
|||
if (l < nct) |
|||
{ |
|||
// Compute the transformation for the l-th column and place the l-th diagonal in VectorS[l].
|
|||
var xnorm = Dnrm2Column(matrixCopy, matrixCopy.RowCount, l, l); |
|||
VectorS[l] = xnorm; |
|||
if (VectorS[l] != 0.0) |
|||
{ |
|||
if (matrixCopy.At(l, l) != 0.0) |
|||
{ |
|||
VectorS[l] = Dsign(VectorS[l], matrixCopy.At(l, l)); |
|||
} |
|||
|
|||
DscalColumn(matrixCopy, matrixCopy.RowCount, l, l, 1.0 / VectorS[l]); |
|||
matrixCopy.At(l, l, (1.0 + matrixCopy.At(l, l))); |
|||
} |
|||
|
|||
VectorS[l] = -VectorS[l]; |
|||
} |
|||
|
|||
for (j = lp1; j < matrixCopy.ColumnCount; j++) |
|||
{ |
|||
if (l < nct) |
|||
{ |
|||
if (VectorS[l] != 0.0) |
|||
{ |
|||
// Apply the transformation.
|
|||
t = -Ddot(matrixCopy, matrixCopy.RowCount, l, j, l) / matrixCopy.At(l, l); |
|||
for (var ii = l; ii < matrixCopy.RowCount; ii++) |
|||
{ |
|||
matrixCopy.At(ii, j, matrixCopy.At(ii, j) + (t * matrixCopy.At(ii, l))); |
|||
} |
|||
} |
|||
} |
|||
|
|||
// Place the l-th row of matrixCopy into e for the
|
|||
// Subsequent calculation of the row transformation.
|
|||
e[j] = matrixCopy.At(l, j); |
|||
} |
|||
|
|||
if (ComputeVectors && l < nct) |
|||
{ |
|||
// Place the transformation in u for subsequent back multiplication.
|
|||
for (i = l; i < matrixCopy.RowCount; i++) |
|||
{ |
|||
MatrixU.At(i, l, matrixCopy.At(i, l)); |
|||
} |
|||
} |
|||
|
|||
if (l >= nrt) |
|||
{ |
|||
continue; |
|||
} |
|||
|
|||
// Compute the l-th row transformation and place the l-th super-diagonal in e(l).
|
|||
var enorm = Dnrm2Vector(e, lp1); |
|||
e[l] = enorm; |
|||
if (e[l] != 0.0) |
|||
{ |
|||
if (e[lp1] != 0.0) |
|||
{ |
|||
e[l] = Dsign(e[l], e[lp1]); |
|||
} |
|||
|
|||
DscalVector(e, lp1, 1.0 / e[l]); |
|||
e[lp1] = 1.0 + e[lp1]; |
|||
} |
|||
|
|||
e[l] = -e[l]; |
|||
if (lp1 < matrixCopy.RowCount && e[l] != 0.0) |
|||
{ |
|||
// Apply the transformation.
|
|||
for (i = lp1; i < matrixCopy.RowCount; i++) |
|||
{ |
|||
work[i] = 0.0; |
|||
} |
|||
|
|||
for (j = lp1; j < matrixCopy.ColumnCount; j++) |
|||
{ |
|||
for (var ii = lp1; ii < matrixCopy.RowCount; ii++) |
|||
{ |
|||
work[ii] += e[j] * matrixCopy.At(ii, j); |
|||
} |
|||
} |
|||
|
|||
for (j = lp1; j < matrixCopy.ColumnCount; j++) |
|||
{ |
|||
var ww = -e[j] / e[lp1]; |
|||
for (var ii = lp1; ii < matrixCopy.RowCount; ii++) |
|||
{ |
|||
matrixCopy.At(ii, j, matrixCopy.At(ii, j) + (ww * work[ii])); |
|||
} |
|||
} |
|||
} |
|||
|
|||
if (ComputeVectors) |
|||
{ |
|||
// Place the transformation in v for subsequent back multiplication.
|
|||
for (i = lp1; i < matrixCopy.ColumnCount; i++) |
|||
{ |
|||
MatrixVT.At(i, l, e[i]); |
|||
} |
|||
} |
|||
} |
|||
|
|||
// Set up the final bidiagonal matrixCopy or order m.
|
|||
var m = Math.Min(matrixCopy.ColumnCount, matrixCopy.RowCount + 1); |
|||
var nctp1 = nct + 1; |
|||
var nrtp1 = nrt + 1; |
|||
if (nct < matrixCopy.ColumnCount) |
|||
{ |
|||
VectorS[nctp1 - 1] = matrixCopy.At((nctp1 - 1), (nctp1 - 1)); |
|||
} |
|||
|
|||
if (matrixCopy.RowCount < m) |
|||
{ |
|||
VectorS[m - 1] = 0.0; |
|||
} |
|||
|
|||
if (nrtp1 < m) |
|||
{ |
|||
e[nrtp1 - 1] = matrixCopy.At((nrtp1 - 1), (m - 1)); |
|||
} |
|||
|
|||
e[m - 1] = 0.0; |
|||
|
|||
// If required, generate u.
|
|||
if (ComputeVectors) |
|||
{ |
|||
for (j = nctp1 - 1; j < ncu; j++) |
|||
{ |
|||
for (i = 0; i < matrixCopy.RowCount; i++) |
|||
{ |
|||
MatrixU.At(i, j, 0.0); |
|||
} |
|||
|
|||
MatrixU.At(j, j, 1.0); |
|||
} |
|||
|
|||
for (l = nct - 1; l >= 0; l--) |
|||
{ |
|||
if (VectorS[l] != 0.0) |
|||
{ |
|||
for (j = l + 1; j < ncu; j++) |
|||
{ |
|||
t = -Ddot(MatrixU, matrixCopy.RowCount, l, j, l) / MatrixU.At(l, l); |
|||
for (var ii = l; ii < matrixCopy.RowCount; ii++) |
|||
{ |
|||
MatrixU.At(ii, j, MatrixU.At(ii, j) + (t * MatrixU.At(ii, l))); |
|||
} |
|||
} |
|||
|
|||
DscalColumn(MatrixU, matrixCopy.RowCount, l, l, -1.0); |
|||
MatrixU.At(l, l, 1.0 + MatrixU.At(l, l)); |
|||
for (i = 0; i < l; i++) |
|||
{ |
|||
MatrixU.At(i, l, 0.0); |
|||
} |
|||
} |
|||
else |
|||
{ |
|||
for (i = 0; i < matrixCopy.RowCount; i++) |
|||
{ |
|||
MatrixU.At(i, l, 0.0); |
|||
} |
|||
|
|||
MatrixU.At(l, l, 1.0); |
|||
} |
|||
} |
|||
} |
|||
|
|||
// If it is required, generate v.
|
|||
if (ComputeVectors) |
|||
{ |
|||
for (l = matrixCopy.ColumnCount - 1; l >= 0; l--) |
|||
{ |
|||
lp1 = l + 1; |
|||
if (l < nrt) |
|||
{ |
|||
if (e[l] != 0.0) |
|||
{ |
|||
for (j = lp1; j < matrixCopy.ColumnCount; j++) |
|||
{ |
|||
t = -Ddot(MatrixVT, matrixCopy.ColumnCount, l, j, lp1) / MatrixVT.At(lp1, l); |
|||
for (var ii = l; ii < matrixCopy.ColumnCount; ii++) |
|||
{ |
|||
MatrixVT.At(ii, j, MatrixVT.At(ii, j) + (t * MatrixVT.At(ii, l))); |
|||
} |
|||
} |
|||
} |
|||
} |
|||
|
|||
for (i = 0; i < matrixCopy.ColumnCount; i++) |
|||
{ |
|||
MatrixVT.At(i, l, 0.0); |
|||
} |
|||
|
|||
MatrixVT.At(l, l, 1.0); |
|||
} |
|||
} |
|||
|
|||
// Transform s and e so that they are double .
|
|||
for (i = 0; i < m; i++) |
|||
{ |
|||
double r; |
|||
if (VectorS[i] != 0.0) |
|||
{ |
|||
t = VectorS[i]; |
|||
r = VectorS[i] / t; |
|||
VectorS[i] = t; |
|||
if (i < m - 1) |
|||
{ |
|||
e[i] = e[i] / r; |
|||
} |
|||
|
|||
if (ComputeVectors) |
|||
{ |
|||
DscalColumn(MatrixU, matrixCopy.RowCount, i, 0, r); |
|||
} |
|||
} |
|||
|
|||
// Exit
|
|||
if (i == m - 1) |
|||
{ |
|||
break; |
|||
} |
|||
|
|||
if (e[i] != 0.0) |
|||
{ |
|||
t = e[i]; |
|||
r = t / e[i]; |
|||
e[i] = t; |
|||
VectorS[i + 1] = VectorS[i + 1] * r; |
|||
if (ComputeVectors) |
|||
{ |
|||
DscalColumn(MatrixVT, matrixCopy.ColumnCount, i + 1, 0, r); |
|||
} |
|||
} |
|||
} |
|||
|
|||
// Main iteration loop for the singular values.
|
|||
var mn = m; |
|||
var iter = 0; |
|||
|
|||
while (m > 0) |
|||
{ |
|||
// Quit if all the singular values have been found. If too many iterations have been performed,
|
|||
// throw exception that Convergence Failed
|
|||
if (iter >= Maxiter) |
|||
{ |
|||
throw new ArgumentException(Resources.ConvergenceFailed); |
|||
} |
|||
|
|||
// This section of the program inspects for negligible elements in the s and e arrays. On
|
|||
// completion the variables kase and l are set as follows.
|
|||
// Kase = 1 if VectorS[m] and e[l-1] are negligible and l < m
|
|||
// Kase = 2 if VectorS[l] is negligible and l < m
|
|||
// Kase = 3 if e[l-1] is negligible, l < m, and VectorS[l, ..., VectorS[m] are not negligible (qr step).
|
|||
// Лase = 4 if e[m-1] is negligible (convergence).
|
|||
double ztest; |
|||
double test; |
|||
for (l = m - 2; l >= 0; l--) |
|||
{ |
|||
test = Math.Abs(VectorS[l]) + Math.Abs(VectorS[l + 1]); |
|||
ztest = test + Math.Abs(e[l]); |
|||
if (ztest.AlmostEqualInDecimalPlaces(test, 15)) |
|||
{ |
|||
e[l] = 0.0; |
|||
break; |
|||
} |
|||
} |
|||
|
|||
int kase; |
|||
if (l == m - 2) |
|||
{ |
|||
kase = 4; |
|||
} |
|||
else |
|||
{ |
|||
int ls; |
|||
for (ls = m - 1; ls > l; ls--) |
|||
{ |
|||
test = 0.0; |
|||
if (ls != m - 1) |
|||
{ |
|||
test = test + Math.Abs(e[ls]); |
|||
} |
|||
|
|||
if (ls != l + 1) |
|||
{ |
|||
test = test + Math.Abs(e[ls - 1]); |
|||
} |
|||
|
|||
ztest = test + Math.Abs(VectorS[ls]); |
|||
if (ztest.AlmostEqualInDecimalPlaces(test, 15)) |
|||
{ |
|||
VectorS[ls] = 0.0; |
|||
break; |
|||
} |
|||
} |
|||
|
|||
if (ls == l) |
|||
{ |
|||
kase = 3; |
|||
} |
|||
else if (ls == m - 1) |
|||
{ |
|||
kase = 1; |
|||
} |
|||
else |
|||
{ |
|||
kase = 2; |
|||
l = ls; |
|||
} |
|||
} |
|||
|
|||
l = l + 1; |
|||
|
|||
// Perform the task indicated by kase.
|
|||
int k; |
|||
double f; |
|||
switch (kase) |
|||
{ |
|||
// Deflate negligible VectorS[m].
|
|||
case 1: |
|||
f = e[m - 2]; |
|||
e[m - 2] = 0.0; |
|||
double t1; |
|||
for (var kk = l; kk < m - 1; kk++) |
|||
{ |
|||
k = m - 2 - kk + l; |
|||
t1 = VectorS[k]; |
|||
Drotg(ref t1, ref f, ref cs, ref sn); |
|||
VectorS[k] = t1; |
|||
if (k != l) |
|||
{ |
|||
f = -sn * e[k - 1]; |
|||
e[k - 1] = cs * e[k - 1]; |
|||
} |
|||
|
|||
if (ComputeVectors) |
|||
{ |
|||
Drot(MatrixVT, matrixCopy.ColumnCount, k, m - 1, cs, sn); |
|||
} |
|||
} |
|||
|
|||
break; |
|||
|
|||
// Split at negligible VectorS[l].
|
|||
case 2: |
|||
f = e[l - 1]; |
|||
e[l - 1] = 0.0; |
|||
for (k = l; k < m; k++) |
|||
{ |
|||
t1 = VectorS[k]; |
|||
Drotg(ref t1, ref f, ref cs, ref sn); |
|||
VectorS[k] = t1; |
|||
f = -sn * e[k]; |
|||
e[k] = cs * e[k]; |
|||
if (ComputeVectors) |
|||
{ |
|||
Drot(MatrixU, matrixCopy.RowCount, k, l - 1, cs, sn); |
|||
} |
|||
} |
|||
|
|||
break; |
|||
|
|||
// Perform one qr step.
|
|||
case 3: |
|||
// Calculate the shift.
|
|||
var scale = 0.0; |
|||
scale = Math.Max(scale, Math.Abs(VectorS[m - 1])); |
|||
scale = Math.Max(scale, Math.Abs(VectorS[m - 2])); |
|||
scale = Math.Max(scale, Math.Abs(e[m - 2])); |
|||
scale = Math.Max(scale, Math.Abs(VectorS[l])); |
|||
scale = Math.Max(scale, Math.Abs(e[l])); |
|||
var sm = VectorS[m - 1] / scale; |
|||
var smm1 = VectorS[m - 2] / scale; |
|||
var emm1 = e[m - 2] / scale; |
|||
var sl = VectorS[l] / scale; |
|||
var el = e[l] / scale; |
|||
var b = (((smm1 + sm) * (smm1 - sm)) + (emm1 * emm1)) / 2.0; |
|||
var c = (sm * emm1) * (sm * emm1); |
|||
var shift = 0.0; |
|||
if (b != 0.0 || c != 0.0) |
|||
{ |
|||
shift = Math.Sqrt((b * b) + c); |
|||
if (b < 0.0) |
|||
{ |
|||
shift = -shift; |
|||
} |
|||
|
|||
shift = c / (b + shift); |
|||
} |
|||
|
|||
f = ((sl + sm) * (sl - sm)) + shift; |
|||
var g = sl * el; |
|||
|
|||
// Chase zeros.
|
|||
for (k = l; k < m - 1; k++) |
|||
{ |
|||
Drotg(ref f, ref g, ref cs, ref sn); |
|||
if (k != l) |
|||
{ |
|||
e[k - 1] = f; |
|||
} |
|||
|
|||
f = (cs * VectorS[k]) + (sn * e[k]); |
|||
e[k] = (cs * e[k]) - (sn * VectorS[k]); |
|||
g = sn * VectorS[k + 1]; |
|||
VectorS[k + 1] = cs * VectorS[k + 1]; |
|||
if (ComputeVectors) |
|||
{ |
|||
Drot(MatrixVT, matrixCopy.ColumnCount, k, k + 1, cs, sn); |
|||
} |
|||
|
|||
Drotg(ref f, ref g, ref cs, ref sn); |
|||
VectorS[k] = f; |
|||
f = (cs * e[k]) + (sn * VectorS[k + 1]); |
|||
VectorS[k + 1] = (-sn * e[k]) + (cs * VectorS[k + 1]); |
|||
g = sn * e[k + 1]; |
|||
e[k + 1] = cs * e[k + 1]; |
|||
if (ComputeVectors && k < matrixCopy.RowCount) |
|||
{ |
|||
Drot(MatrixU, matrixCopy.RowCount, k, k + 1, cs, sn); |
|||
} |
|||
} |
|||
|
|||
e[m - 2] = f; |
|||
iter = iter + 1; |
|||
break; |
|||
|
|||
// Convergence.
|
|||
case 4: |
|||
// Make the singular value positive
|
|||
if (VectorS[l] < 0.0) |
|||
{ |
|||
VectorS[l] = -VectorS[l]; |
|||
if (ComputeVectors) |
|||
{ |
|||
DscalColumn(MatrixVT, matrixCopy.ColumnCount, l, 0, -1.0); |
|||
} |
|||
} |
|||
|
|||
// Order the singular value.
|
|||
while (l != mn - 1) |
|||
{ |
|||
if (VectorS[l] >= VectorS[l + 1]) |
|||
{ |
|||
break; |
|||
} |
|||
|
|||
t = VectorS[l]; |
|||
VectorS[l] = VectorS[l + 1]; |
|||
VectorS[l + 1] = t; |
|||
if (ComputeVectors && l < matrixCopy.ColumnCount) |
|||
{ |
|||
Dswap(MatrixVT, matrixCopy.ColumnCount, l, l + 1); |
|||
} |
|||
|
|||
if (ComputeVectors && l < matrixCopy.RowCount) |
|||
{ |
|||
Dswap(MatrixU, matrixCopy.RowCount, l, l + 1); |
|||
} |
|||
|
|||
l = l + 1; |
|||
} |
|||
|
|||
iter = 0; |
|||
m = m - 1; |
|||
break; |
|||
} |
|||
} |
|||
|
|||
if (ComputeVectors) |
|||
{ |
|||
MatrixVT = MatrixVT.Transpose(); |
|||
} |
|||
|
|||
// Adjust the size of s if rows < columns. We are using ported copy of linpack's svd code and it uses
|
|||
// a singular vector of length mRows+1 when mRows < mColumns. The last element is not used and needs to be removed.
|
|||
// we should port lapack's svd routine to remove this problem.
|
|||
if (matrixCopy.RowCount < matrixCopy.ColumnCount) |
|||
{ |
|||
nm--; |
|||
var tmp = matrixCopy.CreateVector(nm); |
|||
for (i = 0; i < nm; i++) |
|||
{ |
|||
tmp[i] = VectorS[i]; |
|||
} |
|||
|
|||
VectorS = tmp; |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Calculates absolute value of <paramref name="z1"/> multiplied on signum function of <paramref name="z2"/>
|
|||
/// </summary>
|
|||
/// <param name="z1">Double value z1</param>
|
|||
/// <param name="z2">Double value z2</param>
|
|||
/// <returns>Result multiplication of signum function and absolute value</returns>
|
|||
private static double Dsign(double z1, double z2) |
|||
{ |
|||
return Math.Abs(z1) * (z2 / Math.Abs(z2)); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Swap column <paramref name="columnA"/> and <paramref name="columnB"/>
|
|||
/// </summary>
|
|||
/// <param name="a">Source matrix</param>
|
|||
/// <param name="rowCount">The number of rows in <paramref name="a"/></param>
|
|||
/// <param name="columnA">Column A index to swap</param>
|
|||
/// <param name="columnB">Column B index to swap</param>
|
|||
private static void Dswap(Matrix<double> a, int rowCount, int columnA, int columnB) |
|||
{ |
|||
for (var i = 0; i < rowCount; i++) |
|||
{ |
|||
var z = a.At(i, columnA); |
|||
a.At(i, columnA, a.At(i, columnB)); |
|||
a.At(i, columnB, z); |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Scale column <paramref name="column"/> by <paramref name="z"/> starting from row <paramref name="rowStart"/>
|
|||
/// </summary>
|
|||
/// <param name="a">Source matrix</param>
|
|||
/// <param name="rowCount">The number of rows in <paramref name="a"/> </param>
|
|||
/// <param name="column">Column to scale</param>
|
|||
/// <param name="rowStart">Row to scale from</param>
|
|||
/// <param name="z">Scale value</param>
|
|||
private static void DscalColumn(Matrix<double> a, int rowCount, int column, int rowStart, double z) |
|||
{ |
|||
for (var i = rowStart; i < rowCount; i++) |
|||
{ |
|||
a.At(i, column, a.At(i, column) * z); |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Scale vector <paramref name="a"/> by <paramref name="z"/> starting from index <paramref name="start"/>
|
|||
/// </summary>
|
|||
/// <param name="a">Source vector</param>
|
|||
/// <param name="start">Row to scale from</param>
|
|||
/// <param name="z">Scale value</param>
|
|||
private static void DscalVector(double[] a, int start, double z) |
|||
{ |
|||
for (var i = start; i < a.Length; i++) |
|||
{ |
|||
a[i] = a[i] * z; |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Given the Cartesian coordinates (da, db) of a point p, these fucntion return the parameters da, db, c, and s
|
|||
/// associated with the Givens rotation that zeros the y-coordinate of the point.
|
|||
/// </summary>
|
|||
/// <param name="da">Provides the x-coordinate of the point p. On exit contains the parameter r associated with the Givens rotation</param>
|
|||
/// <param name="db">Provides the y-coordinate of the point p. On exit contains the parameter z associated with the Givens rotation</param>
|
|||
/// <param name="c">Contains the parameter c associated with the Givens rotation</param>
|
|||
/// <param name="s">Contains the parameter s associated with the Givens rotation</param>
|
|||
/// <remarks>This is equivalent to the DROTG LAPACK routine.</remarks>
|
|||
private static void Drotg(ref double da, ref double db, ref double c, ref double s) |
|||
{ |
|||
double r, z; |
|||
|
|||
var roe = db; |
|||
var absda = Math.Abs(da); |
|||
var absdb = Math.Abs(db); |
|||
if (absda > absdb) |
|||
{ |
|||
roe = da; |
|||
} |
|||
|
|||
var scale = absda + absdb; |
|||
if (scale == 0.0) |
|||
{ |
|||
c = 1.0; |
|||
s = 0.0; |
|||
r = 0.0; |
|||
z = 0.0; |
|||
} |
|||
else |
|||
{ |
|||
var sda = da / scale; |
|||
var sdb = db / scale; |
|||
r = scale * Math.Sqrt((sda * sda) + (sdb * sdb)); |
|||
if (roe < 0.0) |
|||
{ |
|||
r = -r; |
|||
} |
|||
|
|||
c = da / r; |
|||
s = db / r; |
|||
z = 1.0; |
|||
if (absda > absdb) |
|||
{ |
|||
z = s; |
|||
} |
|||
|
|||
if (absdb >= absda && c != 0.0) |
|||
{ |
|||
z = 1.0 / c; |
|||
} |
|||
} |
|||
|
|||
da = r; |
|||
db = z; |
|||
} |
|||
|
|||
/// <summary>dded
|
|||
/// Calculate Norm 2 of the column <paramref name="column"/> in matrix <paramref name="a"/> starting from row <paramref name="rowStart"/>
|
|||
/// </summary>
|
|||
/// <param name="a">Source matrix</param>
|
|||
/// <param name="rowCount">The number of rows in <paramref name="a"/></param>
|
|||
/// <param name="column">Column index</param>
|
|||
/// <param name="rowStart">Start row index</param>
|
|||
/// <returns>Norm2 (Euclidean norm) of trhe column</returns>
|
|||
private static double Dnrm2Column(Matrix<double> a, int rowCount, int column, int rowStart) |
|||
{ |
|||
double s = 0; |
|||
for (var i = rowStart; i < rowCount; i++) |
|||
{ |
|||
s += a.At(i, column) * a.At(i, column); |
|||
} |
|||
|
|||
return Math.Sqrt(s); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Calculate Norm 2 of the vector <paramref name="a"/> starting from index <paramref name="rowStart"/>
|
|||
/// </summary>
|
|||
/// <param name="a">Source vector</param>
|
|||
/// <param name="rowStart">Start index</param>
|
|||
/// <returns>Norm2 (Euclidean norm) of the vector</returns>
|
|||
private static double Dnrm2Vector(double[] a, int rowStart) |
|||
{ |
|||
double s = 0; |
|||
for (var i = rowStart; i < a.Length; i++) |
|||
{ |
|||
s += a[i] * a[i]; |
|||
} |
|||
|
|||
return Math.Sqrt(s); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Calculate dot product of <paramref name="columnA"/> and <paramref name="columnB"/>
|
|||
/// </summary>
|
|||
/// <param name="a">Source matrix</param>
|
|||
/// <param name="rowCount">The number of rows in <paramref name="a"/></param>
|
|||
/// <param name="columnA">Index of column A</param>
|
|||
/// <param name="columnB">Index of column B</param>
|
|||
/// <param name="rowStart">Starting row index</param>
|
|||
/// <returns>Dot product value</returns>
|
|||
private static double Ddot(Matrix<double> a, int rowCount, int columnA, int columnB, int rowStart) |
|||
{ |
|||
var z = 0.0; |
|||
for (var i = rowStart; i < rowCount; i++) |
|||
{ |
|||
z += a.At(i, columnB) * a.At(i, columnA); |
|||
} |
|||
|
|||
return z; |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Performs rotation of points in the plane. Given two vectors x <paramref name="columnA"/> and y <paramref name="columnB"/>,
|
|||
/// each vector element of these vectors is replaced as follows: x(i) = c*x(i) + s*y(i); y(i) = c*y(i) - s*x(i)
|
|||
/// </summary>
|
|||
/// <param name="a">Source matrix</param>
|
|||
/// <param name="rowCount">The number of rows in <paramref name="a"/></param>
|
|||
/// <param name="columnA">Index of column A</param>
|
|||
/// <param name="columnB">Index of column B</param>
|
|||
/// <param name="c">Scalar "c" value</param>
|
|||
/// <param name="s">Scalar "s" value</param>
|
|||
private static void Drot(Matrix<double> a, int rowCount, int columnA, int columnB, double c, double s) |
|||
{ |
|||
for (var i = 0; i < rowCount; i++) |
|||
{ |
|||
var z = (c * a.At(i, columnA)) + (s * a.At(i, columnB)); |
|||
var tmp = (c * a.At(i, columnB)) - (s * a.At(i, columnA)); |
|||
a.At(i, columnB, tmp); |
|||
a.At(i, columnA, z); |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Solves a system of linear equations, <b>AX = B</b>, with A SVD factorized.
|
|||
/// </summary>
|
|||
/// <param name="input">The right hand side <see cref="Matrix{T}"/>, <b>B</b>.</param>
|
|||
/// <param name="result">The left hand side <see cref="Matrix{T}"/>, <b>X</b>.</param>
|
|||
public override void Solve(Matrix<double> input, Matrix<double> result) |
|||
{ |
|||
// Check for proper arguments.
|
|||
if (input == null) |
|||
{ |
|||
throw new ArgumentNullException("input"); |
|||
} |
|||
|
|||
if (result == null) |
|||
{ |
|||
throw new ArgumentNullException("result"); |
|||
} |
|||
|
|||
if (!ComputeVectors) |
|||
{ |
|||
throw new InvalidOperationException(Resources.SingularVectorsNotComputed); |
|||
} |
|||
|
|||
// The solution X should have the same number of columns as B
|
|||
if (input.ColumnCount != result.ColumnCount) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension); |
|||
} |
|||
|
|||
// The dimension compatibility conditions for X = A\B require the two matrices A and B to have the same number of rows
|
|||
if (MatrixU.RowCount != input.RowCount) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension); |
|||
} |
|||
|
|||
// The solution X row dimension is equal to the column dimension of A
|
|||
if (MatrixVT.ColumnCount != result.RowCount) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension); |
|||
} |
|||
|
|||
var mn = Math.Min(MatrixU.RowCount, MatrixVT.ColumnCount); |
|||
var bn = input.ColumnCount; |
|||
|
|||
var tmp = new double[MatrixVT.ColumnCount]; |
|||
|
|||
for (var k = 0; k < bn; k++) |
|||
{ |
|||
for (var j = 0; j < MatrixVT.ColumnCount; j++) |
|||
{ |
|||
double value = 0; |
|||
if (j < mn) |
|||
{ |
|||
for (var i = 0; i < MatrixU.RowCount; i++) |
|||
{ |
|||
value += MatrixU.At(i, j) * input.At(i, k); |
|||
} |
|||
|
|||
value /= VectorS[j]; |
|||
} |
|||
|
|||
tmp[j] = value; |
|||
} |
|||
|
|||
for (var j = 0; j < MatrixVT.ColumnCount; j++) |
|||
{ |
|||
double value = 0; |
|||
for (var i = 0; i < MatrixVT.ColumnCount; i++) |
|||
{ |
|||
value += MatrixVT.At(i, j) * tmp[i]; |
|||
} |
|||
|
|||
result[j, k] = value; |
|||
} |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Solves a system of linear equations, <b>Ax = b</b>, with A SVD factorized.
|
|||
/// </summary>
|
|||
/// <param name="input">The right hand side vector, <b>b</b>.</param>
|
|||
/// <param name="result">The left hand side <see cref="Matrix{T}"/>, <b>x</b>.</param>
|
|||
public override void Solve(Vector<double> input, Vector<double> result) |
|||
{ |
|||
if (input == null) |
|||
{ |
|||
throw new ArgumentNullException("input"); |
|||
} |
|||
|
|||
if (result == null) |
|||
{ |
|||
throw new ArgumentNullException("result"); |
|||
} |
|||
|
|||
if (!ComputeVectors) |
|||
{ |
|||
throw new InvalidOperationException(Resources.SingularVectorsNotComputed); |
|||
} |
|||
|
|||
// Ax=b where A is an m x n matrix
|
|||
// Check that b is a column vector with m entries
|
|||
if (MatrixU.RowCount != input.Count) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentVectorsSameLength); |
|||
} |
|||
|
|||
// Check that x is a column vector with n entries
|
|||
if (MatrixVT.ColumnCount != result.Count) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentMatrixDimensions); |
|||
} |
|||
|
|||
var mn = Math.Min(MatrixU.RowCount, MatrixVT.ColumnCount); |
|||
var tmp = new double[MatrixVT.ColumnCount]; |
|||
double value; |
|||
for (var j = 0; j < MatrixVT.ColumnCount; j++) |
|||
{ |
|||
value = 0; |
|||
if (j < mn) |
|||
{ |
|||
for (var i = 0; i < MatrixU.RowCount; i++) |
|||
{ |
|||
value += MatrixU.At(i, j) * input[i]; |
|||
} |
|||
|
|||
value /= VectorS[j]; |
|||
} |
|||
|
|||
tmp[j] = value; |
|||
} |
|||
|
|||
for (var j = 0; j < MatrixVT.ColumnCount; j++) |
|||
{ |
|||
value = 0; |
|||
for (int i = 0; i < MatrixVT.ColumnCount; i++) |
|||
{ |
|||
value += MatrixVT.At(i, j) * tmp[i]; |
|||
} |
|||
|
|||
result[j] = value; |
|||
} |
|||
} |
|||
|
|||
#region Simple arithmetic of type T
|
|||
|
|||
/// <summary>
|
|||
/// Multiply two values T*T
|
|||
/// </summary>
|
|||
/// <param name="val1">Left operand value</param>
|
|||
/// <param name="val2">Right operand value</param>
|
|||
/// <returns>Result of multiplication</returns>
|
|||
protected sealed override double MultiplyT(double val1, double val2) |
|||
{ |
|||
return val1 * val2; |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Returns the absolute value of a specified number.
|
|||
/// </summary>
|
|||
/// <param name="val1"> A number whose absolute is to be found</param>
|
|||
/// <returns>Absolute value </returns>
|
|||
protected sealed override double AbsoluteT(double val1) |
|||
{ |
|||
return Math.Abs(val1); |
|||
} |
|||
|
|||
#endregion
|
|||
} |
|||
} |
|||
@ -0,0 +1,118 @@ |
|||
// <copyright file="Svd.cs" company="Math.NET">
|
|||
// Math.NET Numerics, part of the Math.NET Project
|
|||
// http://numerics.mathdotnet.com
|
|||
// http://github.com/mathnet/mathnet-numerics
|
|||
// http://mathnetnumerics.codeplex.com
|
|||
//
|
|||
// Copyright (c) 2009-2010 Math.NET
|
|||
//
|
|||
// Permission is hereby granted, free of charge, to any person
|
|||
// obtaining a copy of this software and associated documentation
|
|||
// files (the "Software"), to deal in the Software without
|
|||
// restriction, including without limitation the rights to use,
|
|||
// copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|||
// copies of the Software, and to permit persons to whom the
|
|||
// Software is furnished to do so, subject to the following
|
|||
// conditions:
|
|||
//
|
|||
// The above copyright notice and this permission notice shall be
|
|||
// included in all copies or substantial portions of the Software.
|
|||
//
|
|||
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
|||
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
|
|||
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
|||
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
|
|||
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
|
|||
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
|||
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
|
|||
// OTHER DEALINGS IN THE SOFTWARE.
|
|||
// </copyright>
|
|||
|
|||
namespace MathNet.Numerics.LinearAlgebra.Double.Factorization |
|||
{ |
|||
using System; |
|||
using System.Linq; |
|||
using Generic; |
|||
using Generic.Factorization; |
|||
using Properties; |
|||
|
|||
/// <summary>
|
|||
/// <para>A class which encapsulates the functionality of the singular value decomposition (SVD).</para>
|
|||
/// <para>Suppose M is an m-by-n matrix whose entries are real numbers.
|
|||
/// Then there exists a factorization of the form M = UΣVT where:
|
|||
/// - U is an m-by-m unitary matrix;
|
|||
/// - Σ is m-by-n diagonal matrix with nonnegative real numbers on the diagonal;
|
|||
/// - VT denotes transpose of V, an n-by-n unitary matrix;
|
|||
/// Such a factorization is called a singular-value decomposition of M. A common convention is to order the diagonal
|
|||
/// entries Σ(i,i) in descending order. In this case, the diagonal matrix Σ is uniquely determined
|
|||
/// by M (though the matrices U and V are not). The diagonal entries of Σ are known as the singular values of M.</para>
|
|||
/// </summary>
|
|||
/// <remarks>
|
|||
/// The computation of the singular value decomposition is done at construction time.
|
|||
/// </remarks>
|
|||
public abstract class Svd : Svd<double> |
|||
{ |
|||
/// <summary>
|
|||
/// Gets the effective numerical matrix rank.
|
|||
/// </summary>
|
|||
/// <value>The number of non-negligible singular values.</value>
|
|||
public override int Rank |
|||
{ |
|||
get |
|||
{ |
|||
return VectorS.Count(t => !Math.Abs(t).AlmostEqual(0.0)); |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Gets the two norm of the <see cref="Matrix{T}"/>.
|
|||
/// </summary>
|
|||
/// <returns>The 2-norm of the <see cref="Matrix{T}"/>.</returns>
|
|||
public override double Norm2 |
|||
{ |
|||
get |
|||
{ |
|||
return Math.Abs(VectorS[0]); |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Gets the condition number <b>max(S) / min(S)</b>
|
|||
/// </summary>
|
|||
/// <returns>The condition number.</returns>
|
|||
public override double ConditionNumber |
|||
{ |
|||
get |
|||
{ |
|||
var tmp = Math.Min(MatrixU.RowCount, MatrixVT.ColumnCount) - 1; |
|||
return Math.Abs(VectorS[0]) / Math.Abs(VectorS[tmp]); |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Gets the determinant of the square matrix for which the SVD was computed.
|
|||
/// </summary>
|
|||
public override double Determinant |
|||
{ |
|||
get |
|||
{ |
|||
if (MatrixU.RowCount != MatrixVT.ColumnCount) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentMatrixSquare); |
|||
} |
|||
|
|||
var det = 1.0; |
|||
foreach (var value in VectorS) |
|||
{ |
|||
det *= value; |
|||
if (Math.Abs(value).AlmostEqual(0.0)) |
|||
{ |
|||
return 0; |
|||
} |
|||
} |
|||
|
|||
return Math.Abs(det); |
|||
} |
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,81 @@ |
|||
// <copyright file="Cholesky.cs" company="Math.NET">
|
|||
// Math.NET Numerics, part of the Math.NET Project
|
|||
// http://numerics.mathdotnet.com
|
|||
// http://github.com/mathnet/mathnet-numerics
|
|||
// http://mathnetnumerics.codeplex.com
|
|||
//
|
|||
// Copyright (c) 2009-2010 Math.NET
|
|||
//
|
|||
// Permission is hereby granted, free of charge, to any person
|
|||
// obtaining a copy of this software and associated documentation
|
|||
// files (the "Software"), to deal in the Software without
|
|||
// restriction, including without limitation the rights to use,
|
|||
// copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|||
// copies of the Software, and to permit persons to whom the
|
|||
// Software is furnished to do so, subject to the following
|
|||
// conditions:
|
|||
//
|
|||
// The above copyright notice and this permission notice shall be
|
|||
// included in all copies or substantial portions of the Software.
|
|||
//
|
|||
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
|||
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
|
|||
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
|||
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
|
|||
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
|
|||
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
|||
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
|
|||
// OTHER DEALINGS IN THE SOFTWARE.
|
|||
// </copyright>
|
|||
|
|||
namespace MathNet.Numerics.LinearAlgebra.Single.Factorization |
|||
{ |
|||
using System; |
|||
using Generic.Factorization; |
|||
|
|||
/// <summary>
|
|||
/// <para>A class which encapsulates the functionality of a Cholesky factorization.</para>
|
|||
/// <para>For a symmetric, positive definite matrix A, the Cholesky factorization
|
|||
/// is an lower triangular matrix L so that A = L*L'.</para>
|
|||
/// </summary>
|
|||
/// <remarks>
|
|||
/// The computation of the Cholesky factorization is done at construction time. If the matrix is not symmetric
|
|||
/// or positive definite, the constructor will throw an exception.
|
|||
/// </remarks>
|
|||
public abstract class Cholesky : Cholesky<float> |
|||
{ |
|||
/// <summary>
|
|||
/// Gets the determinant of the matrix for which the Cholesky matrix was computed.
|
|||
/// </summary>
|
|||
public override float Determinant |
|||
{ |
|||
get |
|||
{ |
|||
var det = 1.0f; |
|||
for (var j = 0; j < CholeskyFactor.RowCount; j++) |
|||
{ |
|||
det *= CholeskyFactor[j, j] * CholeskyFactor[j, j]; |
|||
} |
|||
|
|||
return det; |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Gets the log determinant of the matrix for which the Cholesky matrix was computed.
|
|||
/// </summary>
|
|||
public override float DeterminantLn |
|||
{ |
|||
get |
|||
{ |
|||
var det = 0.0f; |
|||
for (var j = 0; j < CholeskyFactor.RowCount; j++) |
|||
{ |
|||
det += 2.0f * Convert.ToSingle(Math.Log(CholeskyFactor[j, j])); |
|||
} |
|||
|
|||
return det; |
|||
} |
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,115 @@ |
|||
// <copyright file="Evd.cs" company="Math.NET">
|
|||
// Math.NET Numerics, part of the Math.NET Project
|
|||
// http://numerics.mathdotnet.com
|
|||
// http://github.com/mathnet/mathnet-numerics
|
|||
// http://mathnetnumerics.codeplex.com
|
|||
// Copyright (c) 2009-2010 Math.NET
|
|||
// Permission is hereby granted, free of charge, to any person
|
|||
// obtaining a copy of this software and associated documentation
|
|||
// files (the "Software"), to deal in the Software without
|
|||
// restriction, including without limitation the rights to use,
|
|||
// copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|||
// copies of the Software, and to permit persons to whom the
|
|||
// Software is furnished to do so, subject to the following
|
|||
// conditions:
|
|||
// The above copyright notice and this permission notice shall be
|
|||
// included in all copies or substantial portions of the Software.
|
|||
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
|||
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
|
|||
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
|||
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
|
|||
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
|
|||
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
|||
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
|
|||
// OTHER DEALINGS IN THE SOFTWARE.
|
|||
// </copyright>
|
|||
|
|||
namespace MathNet.Numerics.LinearAlgebra.Single.Factorization |
|||
{ |
|||
using System; |
|||
using System.Numerics; |
|||
using Generic.Factorization; |
|||
|
|||
/// <summary>
|
|||
/// Eigenvalues and eigenvectors of a real matrix.
|
|||
/// </summary>
|
|||
/// <remarks>
|
|||
/// If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is
|
|||
/// diagonal and the eigenvector matrix V is orthogonal.
|
|||
/// I.e. A = V*D*V' and V*VT=I.
|
|||
/// If A is not symmetric, then the eigenvalue matrix D is block diagonal
|
|||
/// with the real eigenvalues in 1-by-1 blocks and any complex eigenvalues,
|
|||
/// lambda + i*mu, in 2-by-2 blocks, [lambda, mu; -mu, lambda]. The
|
|||
/// columns of V represent the eigenvectors in the sense that A*V = V*D,
|
|||
/// i.e. A.Multiply(V) equals V.Multiply(D). The matrix V may be badly
|
|||
/// conditioned, or even singular, so the validity of the equation
|
|||
/// A = V*D*Inverse(V) depends upon V.Condition().
|
|||
/// </remarks>
|
|||
public abstract class Evd : Evd<float> |
|||
{ |
|||
/// <summary>
|
|||
/// Gets the absolute value of determinant of the square matrix for which the EVD was computed.
|
|||
/// </summary>
|
|||
public override float Determinant |
|||
{ |
|||
get |
|||
{ |
|||
var det = Complex.One; |
|||
for (var i = 0; i < VectorEv.Count; i++) |
|||
{ |
|||
det *= VectorEv[i]; |
|||
|
|||
if (((Numerics.Complex32)VectorEv[i]).AlmostEqual(Numerics.Complex32.Zero)) |
|||
{ |
|||
return 0; |
|||
} |
|||
} |
|||
|
|||
return Convert.ToSingle(det.Magnitude); |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Gets the effective numerical matrix rank.
|
|||
/// </summary>
|
|||
/// <value>The number of non-negligible singular values.</value>
|
|||
public override int Rank |
|||
{ |
|||
get |
|||
{ |
|||
var rank = 0; |
|||
for (var i = 0; i < VectorEv.Count; i++) |
|||
{ |
|||
if (((Numerics.Complex32)VectorEv[i]).AlmostEqual(Numerics.Complex32.Zero)) |
|||
{ |
|||
continue; |
|||
} |
|||
|
|||
rank++; |
|||
} |
|||
|
|||
return rank; |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Gets a value indicating whether the matrix is full rank or not.
|
|||
/// </summary>
|
|||
/// <value><c>true</c> if the matrix is full rank; otherwise <c>false</c>.</value>
|
|||
public override bool IsFullRank |
|||
{ |
|||
get |
|||
{ |
|||
for (var i = 0; i < VectorEv.Count; i++) |
|||
{ |
|||
if (VectorEv[i].AlmostEqual(Complex.Zero)) |
|||
{ |
|||
return false; |
|||
} |
|||
} |
|||
|
|||
return true; |
|||
} |
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,88 @@ |
|||
// <copyright file="GramSchmidt.cs" company="Math.NET">
|
|||
// Math.NET Numerics, part of the Math.NET Project
|
|||
// http://numerics.mathdotnet.com
|
|||
// http://github.com/mathnet/mathnet-numerics
|
|||
// http://mathnetnumerics.codeplex.com
|
|||
// Copyright (c) 2009-2010 Math.NET
|
|||
// Permission is hereby granted, free of charge, to any person
|
|||
// obtaining a copy of this software and associated documentation
|
|||
// files (the "Software"), to deal in the Software without
|
|||
// restriction, including without limitation the rights to use,
|
|||
// copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|||
// copies of the Software, and to permit persons to whom the
|
|||
// Software is furnished to do so, subject to the following
|
|||
// conditions:
|
|||
// The above copyright notice and this permission notice shall be
|
|||
// included in all copies or substantial portions of the Software.
|
|||
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
|||
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
|
|||
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
|||
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
|
|||
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
|
|||
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
|||
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
|
|||
// OTHER DEALINGS IN THE SOFTWARE.
|
|||
// </copyright>
|
|||
|
|||
namespace MathNet.Numerics.LinearAlgebra.Single.Factorization |
|||
{ |
|||
using System; |
|||
using Generic.Factorization; |
|||
using Properties; |
|||
|
|||
/// <summary>
|
|||
/// <para>A class which encapsulates the functionality of the QR decomposition Modified Gram-Schmidt Orthogonalization.</para>
|
|||
/// <para>Any real square matrix A may be decomposed as A = QR where Q is an orthogonal mxn matrix and R is an nxn upper triangular matrix.</para>
|
|||
/// </summary>
|
|||
/// <remarks>
|
|||
/// The computation of the QR decomposition is done at construction time by modified Gram-Schmidt Orthogonalization.
|
|||
/// </remarks>
|
|||
public abstract class GramSchmidt : GramSchmidt<float> |
|||
{ |
|||
/// <summary>
|
|||
/// Gets the absolute determinant value of the matrix for which the QR matrix was computed.
|
|||
/// </summary>
|
|||
public override float Determinant |
|||
{ |
|||
get |
|||
{ |
|||
if (MatrixR.RowCount != MatrixR.ColumnCount) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentMatrixSquare); |
|||
} |
|||
|
|||
var det = 1.0; |
|||
for (var i = 0; i < MatrixR.ColumnCount; i++) |
|||
{ |
|||
det *= MatrixR.At(i, i); |
|||
if (Math.Abs(MatrixR.At(i, i)).AlmostEqual(0.0f)) |
|||
{ |
|||
return 0; |
|||
} |
|||
} |
|||
|
|||
return Convert.ToSingle(Math.Abs(det)); |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Gets a value indicating whether the matrix is full rank or not.
|
|||
/// </summary>
|
|||
/// <value><c>true</c> if the matrix is full rank; otherwise <c>false</c>.</value>
|
|||
public override bool IsFullRank |
|||
{ |
|||
get |
|||
{ |
|||
for (var i = 0; i < MatrixR.ColumnCount; i++) |
|||
{ |
|||
if (Math.Abs(MatrixR.At(i, i)).AlmostEqual(0.0f)) |
|||
{ |
|||
return false; |
|||
} |
|||
} |
|||
|
|||
return true; |
|||
} |
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,67 @@ |
|||
// <copyright file="LU.cs" company="Math.NET">
|
|||
// Math.NET Numerics, part of the Math.NET Project
|
|||
// http://numerics.mathdotnet.com
|
|||
// http://github.com/mathnet/mathnet-numerics
|
|||
// http://mathnetnumerics.codeplex.com
|
|||
// Copyright (c) 2009-2010 Math.NET
|
|||
// Permission is hereby granted, free of charge, to any person
|
|||
// obtaining a copy of this software and associated documentation
|
|||
// files (the "Software"), to deal in the Software without
|
|||
// restriction, including without limitation the rights to use,
|
|||
// copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|||
// copies of the Software, and to permit persons to whom the
|
|||
// Software is furnished to do so, subject to the following
|
|||
// conditions:
|
|||
// The above copyright notice and this permission notice shall be
|
|||
// included in all copies or substantial portions of the Software.
|
|||
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
|||
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
|
|||
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
|||
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
|
|||
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
|
|||
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
|||
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
|
|||
// OTHER DEALINGS IN THE SOFTWARE.
|
|||
// </copyright>
|
|||
|
|||
namespace MathNet.Numerics.LinearAlgebra.Single.Factorization |
|||
{ |
|||
using Generic.Factorization; |
|||
|
|||
/// <summary>
|
|||
/// <para>A class which encapsulates the functionality of an LU factorization.</para>
|
|||
/// <para>For a matrix A, the LU factorization is a pair of lower triangular matrix L and
|
|||
/// upper triangular matrix U so that A = L*U.</para>
|
|||
/// <para>In the Math.Net implementation we also store a set of pivot elements for increased
|
|||
/// numerical stability. The pivot elements encode a permutation matrix P such that P*A = L*U.</para>
|
|||
/// </summary>
|
|||
/// <remarks>
|
|||
/// The computation of the LU factorization is done at construction time.
|
|||
/// </remarks>
|
|||
public abstract class LU : LU<float> |
|||
{ |
|||
/// <summary>
|
|||
/// Gets the determinant of the matrix for which the LU factorization was computed.
|
|||
/// </summary>
|
|||
public override float Determinant |
|||
{ |
|||
get |
|||
{ |
|||
var det = 1.0f; |
|||
for (var j = 0; j < Factors.RowCount; j++) |
|||
{ |
|||
if (Pivots[j] != j) |
|||
{ |
|||
det *= -Factors.At(j, j); |
|||
} |
|||
else |
|||
{ |
|||
det *= Factors.At(j, j); |
|||
} |
|||
} |
|||
|
|||
return det; |
|||
} |
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,90 @@ |
|||
// <copyright file="QR.cs" company="Math.NET">
|
|||
// Math.NET Numerics, part of the Math.NET Project
|
|||
// http://numerics.mathdotnet.com
|
|||
// http://github.com/mathnet/mathnet-numerics
|
|||
// http://mathnetnumerics.codeplex.com
|
|||
// Copyright (c) 2009-2010 Math.NET
|
|||
// Permission is hereby granted, free of charge, to any person
|
|||
// obtaining a copy of this software and associated documentation
|
|||
// files (the "Software"), to deal in the Software without
|
|||
// restriction, including without limitation the rights to use,
|
|||
// copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|||
// copies of the Software, and to permit persons to whom the
|
|||
// Software is furnished to do so, subject to the following
|
|||
// conditions:
|
|||
// The above copyright notice and this permission notice shall be
|
|||
// included in all copies or substantial portions of the Software.
|
|||
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
|||
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
|
|||
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
|||
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
|
|||
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
|
|||
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
|||
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
|
|||
// OTHER DEALINGS IN THE SOFTWARE.
|
|||
// </copyright>
|
|||
|
|||
namespace MathNet.Numerics.LinearAlgebra.Single.Factorization |
|||
{ |
|||
using System; |
|||
using Generic.Factorization; |
|||
using Properties; |
|||
|
|||
/// <summary>
|
|||
/// <para>A class which encapsulates the functionality of the QR decomposition.</para>
|
|||
/// <para>Any real square matrix A (m x n) may be decomposed as A = QR where Q is an orthogonal matrix (m x m)
|
|||
/// (its columns are orthogonal unit vectors meaning QTQ = I) and R (m x n) is an upper triangular matrix
|
|||
/// (also called right triangular matrix).</para>
|
|||
/// </summary>
|
|||
/// <remarks>
|
|||
/// The computation of the QR decomposition is done at construction time by Householder transformation.
|
|||
/// </remarks>
|
|||
public abstract class QR : QR<float> |
|||
{ |
|||
/// <summary>
|
|||
/// Gets the absolute determinant value of the matrix for which the QR matrix was computed.
|
|||
/// </summary>
|
|||
public override float Determinant |
|||
{ |
|||
get |
|||
{ |
|||
if (MatrixR.RowCount != MatrixR.ColumnCount) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentMatrixSquare); |
|||
} |
|||
|
|||
var det = 1.0; |
|||
for (var i = 0; i < MatrixR.ColumnCount; i++) |
|||
{ |
|||
det *= MatrixR.At(i, i); |
|||
if (Math.Abs(MatrixR.At(i, i)).AlmostEqual(0.0f)) |
|||
{ |
|||
return 0; |
|||
} |
|||
} |
|||
|
|||
return Convert.ToSingle(Math.Abs(det)); |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Gets a value indicating whether the matrix is full rank or not.
|
|||
/// </summary>
|
|||
/// <value><c>true</c> if the matrix is full rank; otherwise <c>false</c>.</value>
|
|||
public override bool IsFullRank |
|||
{ |
|||
get |
|||
{ |
|||
for (var i = 0; i < MatrixR.ColumnCount; i++) |
|||
{ |
|||
if (Math.Abs(MatrixR.At(i, i)).AlmostEqual(0.0f)) |
|||
{ |
|||
return false; |
|||
} |
|||
} |
|||
|
|||
return true; |
|||
} |
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,118 @@ |
|||
// <copyright file="Svd.cs" company="Math.NET">
|
|||
// Math.NET Numerics, part of the Math.NET Project
|
|||
// http://numerics.mathdotnet.com
|
|||
// http://github.com/mathnet/mathnet-numerics
|
|||
// http://mathnetnumerics.codeplex.com
|
|||
//
|
|||
// Copyright (c) 2009-2010 Math.NET
|
|||
//
|
|||
// Permission is hereby granted, free of charge, to any person
|
|||
// obtaining a copy of this software and associated documentation
|
|||
// files (the "Software"), to deal in the Software without
|
|||
// restriction, including without limitation the rights to use,
|
|||
// copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|||
// copies of the Software, and to permit persons to whom the
|
|||
// Software is furnished to do so, subject to the following
|
|||
// conditions:
|
|||
//
|
|||
// The above copyright notice and this permission notice shall be
|
|||
// included in all copies or substantial portions of the Software.
|
|||
//
|
|||
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
|||
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
|
|||
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
|||
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
|
|||
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
|
|||
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
|||
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
|
|||
// OTHER DEALINGS IN THE SOFTWARE.
|
|||
// </copyright>
|
|||
|
|||
namespace MathNet.Numerics.LinearAlgebra.Single.Factorization |
|||
{ |
|||
using System; |
|||
using System.Linq; |
|||
using Generic; |
|||
using Generic.Factorization; |
|||
using Properties; |
|||
|
|||
/// <summary>
|
|||
/// <para>A class which encapsulates the functionality of the singular value decomposition (SVD).</para>
|
|||
/// <para>Suppose M is an m-by-n matrix whose entries are real numbers.
|
|||
/// Then there exists a factorization of the form M = UΣVT where:
|
|||
/// - U is an m-by-m unitary matrix;
|
|||
/// - Σ is m-by-n diagonal matrix with nonnegative real numbers on the diagonal;
|
|||
/// - VT denotes transpose of V, an n-by-n unitary matrix;
|
|||
/// Such a factorization is called a singular-value decomposition of M. A common convention is to order the diagonal
|
|||
/// entries Σ(i,i) in descending order. In this case, the diagonal matrix Σ is uniquely determined
|
|||
/// by M (though the matrices U and V are not). The diagonal entries of Σ are known as the singular values of M.</para>
|
|||
/// </summary>
|
|||
/// <remarks>
|
|||
/// The computation of the singular value decomposition is done at construction time.
|
|||
/// </remarks>
|
|||
public abstract class Svd : Svd<float> |
|||
{ |
|||
/// <summary>
|
|||
/// Gets the effective numerical matrix rank.
|
|||
/// </summary>
|
|||
/// <value>The number of non-negligible singular values.</value>
|
|||
public override int Rank |
|||
{ |
|||
get |
|||
{ |
|||
return VectorS.Count(t => !Math.Abs(t).AlmostEqual(0.0f)); |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Gets the two norm of the <see cref="Matrix{T}"/>.
|
|||
/// </summary>
|
|||
/// <returns>The 2-norm of the <see cref="Matrix{T}"/>.</returns>
|
|||
public override float Norm2 |
|||
{ |
|||
get |
|||
{ |
|||
return Math.Abs(VectorS[0]); |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Gets the condition number <b>max(S) / min(S)</b>
|
|||
/// </summary>
|
|||
/// <returns>The condition number.</returns>
|
|||
public override float ConditionNumber |
|||
{ |
|||
get |
|||
{ |
|||
var tmp = Math.Min(MatrixU.RowCount, MatrixVT.ColumnCount) - 1; |
|||
return Math.Abs(VectorS[0]) / Math.Abs(VectorS[tmp]); |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Gets the determinant of the square matrix for which the SVD was computed.
|
|||
/// </summary>
|
|||
public override float Determinant |
|||
{ |
|||
get |
|||
{ |
|||
if (MatrixU.RowCount != MatrixVT.ColumnCount) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentMatrixSquare); |
|||
} |
|||
|
|||
var det = 1.0; |
|||
foreach (var value in VectorS) |
|||
{ |
|||
det *= value; |
|||
if (Math.Abs(value).AlmostEqual(0.0f)) |
|||
{ |
|||
return 0; |
|||
} |
|||
} |
|||
|
|||
return Convert.ToSingle(Math.Abs(det)); |
|||
} |
|||
} |
|||
} |
|||
} |
|||
Loading…
Reference in new issue