Browse Source

LA/Precision: Matrix norms double-valued (like vectors); generic AlmostEqual support

pull/163/merge
Christoph Ruegg 13 years ago
parent
commit
59b1005d7f
  1. 2
      src/Examples/LinearAlgebra/Factorization/Svd.cs
  2. 20
      src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs
  3. 21
      src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs
  4. 2
      src/Numerics/LinearAlgebra/Complex/Factorization/Svd.cs
  5. 60
      src/Numerics/LinearAlgebra/Complex/Matrix.cs
  6. 52
      src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs
  7. 20
      src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs
  8. 21
      src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs
  9. 2
      src/Numerics/LinearAlgebra/Complex32/Factorization/Svd.cs
  10. 61
      src/Numerics/LinearAlgebra/Complex32/Matrix.cs
  11. 56
      src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs
  12. 14
      src/Numerics/LinearAlgebra/Double/DenseMatrix.cs
  13. 15
      src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs
  14. 2
      src/Numerics/LinearAlgebra/Double/Factorization/Svd.cs
  15. 56
      src/Numerics/LinearAlgebra/Double/Matrix.cs
  16. 53
      src/Numerics/LinearAlgebra/Double/SparseMatrix.cs
  17. 2
      src/Numerics/LinearAlgebra/Factorization/Svd.cs
  18. 21
      src/Numerics/LinearAlgebra/Matrix.Arithmetic.cs
  19. 21
      src/Numerics/LinearAlgebra/Matrix.cs
  20. 18
      src/Numerics/LinearAlgebra/Single/DenseMatrix.cs
  21. 21
      src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs
  22. 2
      src/Numerics/LinearAlgebra/Single/Factorization/Svd.cs
  23. 60
      src/Numerics/LinearAlgebra/Single/Matrix.cs
  24. 57
      src/Numerics/LinearAlgebra/Single/SparseMatrix.cs
  25. 138
      src/Numerics/Precision.Equality.cs
  26. 25
      src/UnitTests/LinearAlgebraTests/Complex32/MatrixTests.cs

2
src/Examples/LinearAlgebra/Factorization/Svd.cs

@ -135,7 +135,7 @@ namespace Examples.LinearAlgebra.FactorizationExamples
// 10. 2-norm of the matrix // 10. 2-norm of the matrix
Console.WriteLine(@"10. 2-norm of the matrix"); Console.WriteLine(@"10. 2-norm of the matrix");
Console.WriteLine(svd.Norm2); Console.WriteLine(svd.L2Norm);
Console.WriteLine(); Console.WriteLine();
// 11. Rank of the matrix // 11. Rank of the matrix

20
src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs

@ -428,23 +428,23 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// <summary>Calculates the L1 norm.</summary> /// <summary>Calculates the L1 norm.</summary>
/// <returns>The L1 norm of the matrix.</returns> /// <returns>The L1 norm of the matrix.</returns>
public override Complex L1Norm() public override double L1Norm()
{ {
return Control.LinearAlgebraProvider.MatrixNorm(Norm.OneNorm, _rowCount, _columnCount, _values); return Control.LinearAlgebraProvider.MatrixNorm(Norm.OneNorm, _rowCount, _columnCount, _values).Real;
} }
/// <summary>Calculates the Frobenius norm of this matrix.</summary> /// <summary>Calculates the infinity norm of this matrix.</summary>
/// <returns>The Frobenius norm of this matrix.</returns> /// <returns>The infinity norm of this matrix.</returns>
public override Complex FrobeniusNorm() public override double InfinityNorm()
{ {
return Control.LinearAlgebraProvider.MatrixNorm(Norm.FrobeniusNorm, _rowCount, _columnCount, _values); return Control.LinearAlgebraProvider.MatrixNorm(Norm.InfinityNorm, _rowCount, _columnCount, _values).Real;
} }
/// <summary>Calculates the infinity norm of this matrix.</summary> /// <summary>Calculates the Frobenius norm of this matrix.</summary>
/// <returns>The infinity norm of this matrix.</returns> /// <returns>The Frobenius norm of this matrix.</returns>
public override Complex InfinityNorm() public override double FrobeniusNorm()
{ {
return Control.LinearAlgebraProvider.MatrixNorm(Norm.InfinityNorm, _rowCount, _columnCount, _values); return Control.LinearAlgebraProvider.MatrixNorm(Norm.FrobeniusNorm, _rowCount, _columnCount, _values).Real;
} }
/// <summary> /// <summary>

21
src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs

@ -654,33 +654,32 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// <summary>Calculates the L1 norm.</summary> /// <summary>Calculates the L1 norm.</summary>
/// <returns>The L1 norm of the matrix.</returns> /// <returns>The L1 norm of the matrix.</returns>
public override Complex L1Norm() public override double L1Norm()
{ {
return _data.Aggregate(double.NegativeInfinity, (current, t) => Math.Max(current, t.Magnitude)); return _data.Aggregate(double.NegativeInfinity, (current, t) => Math.Max(current, t.Magnitude));
} }
/// <summary>Calculates the L2 norm.</summary> /// <summary>Calculates the L2 norm.</summary>
/// <returns>The L2 norm of the matrix.</returns> /// <returns>The L2 norm of the matrix.</returns>
public override Complex L2Norm() public override double L2Norm()
{ {
return _data.Aggregate(double.NegativeInfinity, (current, t) => Math.Max(current, t.Magnitude)); return _data.Aggregate(double.NegativeInfinity, (current, t) => Math.Max(current, t.Magnitude));
} }
/// <summary>Calculates the Frobenius norm of this matrix.</summary>
/// <returns>The Frobenius norm of this matrix.</returns>
public override Complex FrobeniusNorm()
{
var norm = _data.Sum(t => t.Magnitude * t.Magnitude);
return Math.Sqrt(norm);
}
/// <summary>Calculates the infinity norm of this matrix.</summary> /// <summary>Calculates the infinity norm of this matrix.</summary>
/// <returns>The infinity norm of this matrix.</returns> /// <returns>The infinity norm of this matrix.</returns>
public override Complex InfinityNorm() public override double InfinityNorm()
{ {
return L1Norm(); return L1Norm();
} }
/// <summary>Calculates the Frobenius norm of this matrix.</summary>
/// <returns>The Frobenius norm of this matrix.</returns>
public override double FrobeniusNorm()
{
return Math.Sqrt(_data.Sum(t => t.Magnitude * t.Magnitude));
}
/// <summary>Calculates the condition number of this matrix.</summary> /// <summary>Calculates the condition number of this matrix.</summary>
/// <returns>The condition number of the matrix.</returns> /// <returns>The condition number of the matrix.</returns>
public override Complex ConditionNumber() public override Complex ConditionNumber()

2
src/Numerics/LinearAlgebra/Complex/Factorization/Svd.cs

@ -79,7 +79,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
/// Gets the two norm of the <see cref="Matrix{T}"/>. /// Gets the two norm of the <see cref="Matrix{T}"/>.
/// </summary> /// </summary>
/// <returns>The 2-norm of the <see cref="Matrix{T}"/>.</returns> /// <returns>The 2-norm of the <see cref="Matrix{T}"/>.</returns>
public override Complex Norm2 public override double L2Norm
{ {
get get
{ {

60
src/Numerics/LinearAlgebra/Complex/Matrix.cs

@ -59,76 +59,68 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// <summary>Calculates the L1 norm.</summary> /// <summary>Calculates the L1 norm.</summary>
/// <returns>The L1 norm of the matrix.</returns> /// <returns>The L1 norm of the matrix.</returns>
public override Complex L1Norm() public override double L1Norm()
{ {
var norm = 0.0; var norm = 0d;
for (var j = 0; j < ColumnCount; j++) for (var j = 0; j < ColumnCount; j++)
{ {
var s = 0.0; var s = 0d;
for (var i = 0; i < RowCount; i++) for (var i = 0; i < RowCount; i++)
{ {
s += At(i, j).Magnitude; s += At(i, j).Magnitude;
} }
norm = Math.Max(norm, s); norm = Math.Max(norm, s);
} }
return norm; return norm;
} }
/// <summary> /// <summary>Calculates the infinity norm of this matrix.</summary>
/// Returns the conjugate transpose of this matrix. /// <returns>The infinity norm of this matrix.</returns>
/// </summary> public override double InfinityNorm()
/// <returns>The conjugate transpose of this matrix.</returns>
public override Matrix<Complex> ConjugateTranspose()
{ {
var ret = CreateMatrix(ColumnCount, RowCount); var norm = 0d;
for (var j = 0; j < ColumnCount; j++) for (var i = 0; i < RowCount; i++)
{ {
for (var i = 0; i < RowCount; i++) var s = 0d;
for (var j = 0; j < ColumnCount; j++)
{ {
ret.At(j, i, At(i, j).Conjugate()); s += At(i, j).Magnitude;
} }
norm = Math.Max(norm, s);
} }
return norm;
return ret;
} }
/// <summary>Calculates the Frobenius norm of this matrix.</summary> /// <summary>Calculates the Frobenius norm of this matrix.</summary>
/// <returns>The Frobenius norm of this matrix.</returns> /// <returns>The Frobenius norm of this matrix.</returns>
public override Complex FrobeniusNorm() public override double FrobeniusNorm()
{ {
var transpose = ConjugateTranspose(); var transpose = ConjugateTranspose();
var aat = this * transpose; var aat = this * transpose;
var norm = 0d;
var norm = 0.0;
for (var i = 0; i < RowCount; i++) for (var i = 0; i < RowCount; i++)
{ {
norm += aat.At(i, i).Magnitude; norm += aat.At(i, i).Magnitude;
} }
return Math.Sqrt(norm);
norm = Math.Sqrt(norm);
return norm;
} }
/// <summary>Calculates the infinity norm of this matrix.</summary> /// <summary>
/// <returns>The infinity norm of this matrix.</returns> /// Returns the conjugate transpose of this matrix.
public override Complex InfinityNorm() /// </summary>
/// <returns>The conjugate transpose of this matrix.</returns>
public override Matrix<Complex> ConjugateTranspose()
{ {
var norm = 0.0; var ret = CreateMatrix(ColumnCount, RowCount);
for (var i = 0; i < RowCount; i++) for (var j = 0; j < ColumnCount; j++)
{ {
var s = 0.0; for (var i = 0; i < RowCount; i++)
for (var j = 0; j < ColumnCount; j++)
{ {
s += At(i, j).Magnitude; ret.At(j, i, At(i, j).Conjugate());
} }
norm = Math.Max(norm, s);
} }
return norm; return ret;
} }
/// <summary> /// <summary>

52
src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs

@ -674,64 +674,58 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
return new SparseMatrix(ret); return new SparseMatrix(ret);
} }
/// <summary>Calculates the Frobenius norm of this matrix.</summary> /// <summary>Calculates the infinity norm of this matrix.</summary>
/// <returns>The Frobenius norm of this matrix.</returns> /// <returns>The infinity norm of this matrix.</returns>
public override Complex FrobeniusNorm() public override double InfinityNorm()
{ {
var aat = (SparseCompressedRowMatrixStorage<Complex>) (this*ConjugateTranspose()).Storage; var rowPointers = _storage.RowPointers;
var values = _storage.Values;
var norm = 0d; var norm = 0d;
for (var i = 0; i < RowCount; i++)
for (var i = 0; i < aat.RowCount; i++)
{ {
var startIndex = aat.RowPointers[i]; var startIndex = rowPointers[i];
var endIndex = aat.RowPointers[i + 1]; var endIndex = rowPointers[i + 1];
if (startIndex == endIndex) if (startIndex == endIndex)
{ {
continue; continue;
} }
var s = 0d;
for (var j = startIndex; j < endIndex; j++) for (var j = startIndex; j < endIndex; j++)
{ {
if (i == aat.ColumnIndices[j]) s += values[j].Magnitude;
{
norm += aat.Values[j].Magnitude;
}
} }
norm = Math.Max(norm, s);
} }
norm = Math.Sqrt(norm);
return norm; return norm;
} }
/// <summary>Calculates the infinity norm of this matrix.</summary> /// <summary>Calculates the Frobenius norm of this matrix.</summary>
/// <returns>The infinity norm of this matrix.</returns> /// <returns>The Frobenius norm of this matrix.</returns>
public override Complex InfinityNorm() public override double FrobeniusNorm()
{ {
var rowPointers = _storage.RowPointers; var aat = (SparseCompressedRowMatrixStorage<Complex>) (this*ConjugateTranspose()).Storage;
var values = _storage.Values;
var norm = 0d; var norm = 0d;
for (var i = 0; i < RowCount; i++) for (var i = 0; i < aat.RowCount; i++)
{ {
var startIndex = rowPointers[i]; var startIndex = aat.RowPointers[i];
var endIndex = rowPointers[i + 1]; var endIndex = aat.RowPointers[i + 1];
if (startIndex == endIndex) if (startIndex == endIndex)
{ {
continue; continue;
} }
var s = 0.0;
for (var j = startIndex; j < endIndex; j++) for (var j = startIndex; j < endIndex; j++)
{ {
s += values[j].Magnitude; if (i == aat.ColumnIndices[j])
{
norm += aat.Values[j].Magnitude;
}
} }
norm = Math.Max(norm, s);
} }
return Math.Sqrt(norm);
return norm;
} }
/// <summary> /// <summary>

20
src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs

@ -423,23 +423,23 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// <summary>Calculates the L1 norm.</summary> /// <summary>Calculates the L1 norm.</summary>
/// <returns>The L1 norm of the matrix.</returns> /// <returns>The L1 norm of the matrix.</returns>
public override Complex32 L1Norm() public override double L1Norm()
{ {
return Control.LinearAlgebraProvider.MatrixNorm(Norm.OneNorm, _rowCount, _columnCount, _values); return Control.LinearAlgebraProvider.MatrixNorm(Norm.OneNorm, _rowCount, _columnCount, _values).Real;
} }
/// <summary>Calculates the Frobenius norm of this matrix.</summary> /// <summary>Calculates the infinity norm of this matrix.</summary>
/// <returns>The Frobenius norm of this matrix.</returns> /// <returns>The infinity norm of this matrix.</returns>
public override Complex32 FrobeniusNorm() public override double InfinityNorm()
{ {
return Control.LinearAlgebraProvider.MatrixNorm(Norm.FrobeniusNorm, _rowCount, _columnCount, _values); return Control.LinearAlgebraProvider.MatrixNorm(Norm.InfinityNorm, _rowCount, _columnCount, _values).Real;
} }
/// <summary>Calculates the infinity norm of this matrix.</summary> /// <summary>Calculates the Frobenius norm of this matrix.</summary>
/// <returns>The infinity norm of this matrix.</returns> /// <returns>The Frobenius norm of this matrix.</returns>
public override Complex32 InfinityNorm() public override double FrobeniusNorm()
{ {
return Control.LinearAlgebraProvider.MatrixNorm(Norm.InfinityNorm, _rowCount, _columnCount, _values); return Control.LinearAlgebraProvider.MatrixNorm(Norm.FrobeniusNorm, _rowCount, _columnCount, _values).Real;
} }
/// <summary> /// <summary>

21
src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs

@ -649,33 +649,32 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// <summary>Calculates the L1 norm.</summary> /// <summary>Calculates the L1 norm.</summary>
/// <returns>The L1 norm of the matrix.</returns> /// <returns>The L1 norm of the matrix.</returns>
public override Complex32 L1Norm() public override double L1Norm()
{ {
return _data.Aggregate(float.NegativeInfinity, (current, t) => Math.Max(current, t.Magnitude)); return _data.Aggregate(float.NegativeInfinity, (current, t) => Math.Max(current, t.Magnitude));
} }
/// <summary>Calculates the L2 norm.</summary> /// <summary>Calculates the L2 norm.</summary>
/// <returns>The L2 norm of the matrix.</returns> /// <returns>The L2 norm of the matrix.</returns>
public override Complex32 L2Norm() public override double L2Norm()
{ {
return _data.Aggregate(float.NegativeInfinity, (current, t) => Math.Max(current, t.Magnitude)); return _data.Aggregate(float.NegativeInfinity, (current, t) => Math.Max(current, t.Magnitude));
} }
/// <summary>Calculates the Frobenius norm of this matrix.</summary>
/// <returns>The Frobenius norm of this matrix.</returns>
public override Complex32 FrobeniusNorm()
{
var norm = _data.Sum(t => t.Magnitude * t.Magnitude);
return Convert.ToSingle(Math.Sqrt(norm));
}
/// <summary>Calculates the infinity norm of this matrix.</summary> /// <summary>Calculates the infinity norm of this matrix.</summary>
/// <returns>The infinity norm of this matrix.</returns> /// <returns>The infinity norm of this matrix.</returns>
public override Complex32 InfinityNorm() public override double InfinityNorm()
{ {
return L1Norm(); return L1Norm();
} }
/// <summary>Calculates the Frobenius norm of this matrix.</summary>
/// <returns>The Frobenius norm of this matrix.</returns>
public override double FrobeniusNorm()
{
return Math.Sqrt(_data.Sum(t => t.Magnitude * t.Magnitude));
}
/// <summary>Calculates the condition number of this matrix.</summary> /// <summary>Calculates the condition number of this matrix.</summary>
/// <returns>The condition number of the matrix.</returns> /// <returns>The condition number of the matrix.</returns>
public override Complex32 ConditionNumber() public override Complex32 ConditionNumber()

2
src/Numerics/LinearAlgebra/Complex32/Factorization/Svd.cs

@ -74,7 +74,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
/// Gets the two norm of the <see cref="Matrix{T}"/>. /// Gets the two norm of the <see cref="Matrix{T}"/>.
/// </summary> /// </summary>
/// <returns>The 2-norm of the <see cref="Matrix{T}"/>.</returns> /// <returns>The 2-norm of the <see cref="Matrix{T}"/>.</returns>
public override Complex32 Norm2 public override double L2Norm
{ {
get get
{ {

61
src/Numerics/LinearAlgebra/Complex32/Matrix.cs

@ -54,76 +54,67 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// <summary>Calculates the L1 norm.</summary> /// <summary>Calculates the L1 norm.</summary>
/// <returns>The L1 norm of the matrix.</returns> /// <returns>The L1 norm of the matrix.</returns>
public override Complex32 L1Norm() public override double L1Norm()
{ {
var norm = 0.0f; var norm = 0d;
for (var j = 0; j < ColumnCount; j++) for (var j = 0; j < ColumnCount; j++)
{ {
var s = 0.0f; var s = 0d;
for (var i = 0; i < RowCount; i++) for (var i = 0; i < RowCount; i++)
{ {
s += At(i, j).Magnitude; s += At(i, j).Magnitude;
} }
norm = Math.Max(norm, s); norm = Math.Max(norm, s);
} }
return norm; return norm;
} }
/// <summary> /// <summary>Calculates the infinity norm of this matrix.</summary>
/// Returns the conjugate transpose of this matrix. /// <returns>The infinity norm of this matrix.</returns>
/// </summary> public override double InfinityNorm()
/// <returns>The conjugate transpose of this matrix.</returns>
public override Matrix<Complex32> ConjugateTranspose()
{ {
var ret = CreateMatrix(ColumnCount, RowCount); var norm = 0d;
for (var j = 0; j < ColumnCount; j++) for (var i = 0; i < RowCount; i++)
{ {
for (var i = 0; i < RowCount; i++) var s = 0d;
for (var j = 0; j < ColumnCount; j++)
{ {
ret.At(j, i, At(i, j).Conjugate()); s += At(i, j).Magnitude;
} }
norm = Math.Max(norm, s);
} }
return norm;
return ret;
} }
/// <summary>Calculates the Frobenius norm of this matrix.</summary> /// <summary>Calculates the Frobenius norm of this matrix.</summary>
/// <returns>The Frobenius norm of this matrix.</returns> /// <returns>The Frobenius norm of this matrix.</returns>
public override Complex32 FrobeniusNorm() public override double FrobeniusNorm()
{ {
var transpose = ConjugateTranspose(); var transpose = ConjugateTranspose();
var aat = this * transpose; var aat = this * transpose;
var norm = 0d;
var norm = 0.0f;
for (var i = 0; i < RowCount; i++) for (var i = 0; i < RowCount; i++)
{ {
norm += aat.At(i, i).Magnitude; norm += aat.At(i, i).Magnitude;
} }
return Math.Sqrt(norm);
norm = Convert.ToSingle(Math.Sqrt(norm));
return norm;
} }
/// <summary>Calculates the infinity norm of this matrix.</summary> /// <summary>
/// <returns>The infinity norm of this matrix.</returns> /// Returns the conjugate transpose of this matrix.
public override Complex32 InfinityNorm() /// </summary>
/// <returns>The conjugate transpose of this matrix.</returns>
public override Matrix<Complex32> ConjugateTranspose()
{ {
var norm = 0.0f; var ret = CreateMatrix(ColumnCount, RowCount);
for (var i = 0; i < RowCount; i++) for (var j = 0; j < ColumnCount; j++)
{ {
var s = 0.0f; for (var i = 0; i < RowCount; i++)
for (var j = 0; j < ColumnCount; j++)
{ {
s += At(i, j).Magnitude; ret.At(j, i, At(i, j).Conjugate());
} }
norm = Math.Max(norm, s);
} }
return ret;
return norm;
} }
/// <summary> /// <summary>

56
src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs

@ -669,64 +669,58 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
return new SparseMatrix(ret); return new SparseMatrix(ret);
} }
/// <summary>Calculates the Frobenius norm of this matrix.</summary> /// <summary>Calculates the infinity norm of this matrix.</summary>
/// <returns>The Frobenius norm of this matrix.</returns> /// <returns>The infinity norm of this matrix.</returns>
public override Complex32 FrobeniusNorm() public override double InfinityNorm()
{ {
var aat = (SparseCompressedRowMatrixStorage<Complex32>) (this*ConjugateTranspose()).Storage; var rowPointers = _storage.RowPointers;
var norm = 0f; var values = _storage.Values;
var norm = 0d;
for (var i = 0; i < aat.RowCount; i++) for (var i = 0; i < RowCount; i++)
{ {
var startIndex = aat.RowPointers[i]; var startIndex = rowPointers[i];
var endIndex = aat.RowPointers[i + 1]; var endIndex = rowPointers[i + 1];
if (startIndex == endIndex) if (startIndex == endIndex)
{ {
continue; continue;
} }
var s = 0d;
for (var j = startIndex; j < endIndex; j++) for (var j = startIndex; j < endIndex; j++)
{ {
if (i == aat.ColumnIndices[j]) s += values[j].Magnitude;
{
norm += aat.Values[j].Magnitude;
}
} }
norm = Math.Max(norm, s);
} }
norm = Convert.ToSingle(Math.Sqrt(norm));
return norm; return norm;
} }
/// <summary>Calculates the infinity norm of this matrix.</summary> /// <summary>Calculates the Frobenius norm of this matrix.</summary>
/// <returns>The infinity norm of this matrix.</returns> /// <returns>The Frobenius norm of this matrix.</returns>
public override Complex32 InfinityNorm() public override double FrobeniusNorm()
{ {
var rowPointers = _storage.RowPointers; var aat = (SparseCompressedRowMatrixStorage<Complex32>) (this*ConjugateTranspose()).Storage;
var values = _storage.Values; var norm = 0d;
for (var i = 0; i < aat.RowCount; i++)
var norm = 0f;
for (var i = 0; i < RowCount; i++)
{ {
var startIndex = rowPointers[i]; var startIndex = aat.RowPointers[i];
var endIndex = rowPointers[i + 1]; var endIndex = aat.RowPointers[i + 1];
if (startIndex == endIndex) if (startIndex == endIndex)
{ {
continue; continue;
} }
var s = 0.0f;
for (var j = startIndex; j < endIndex; j++) for (var j = startIndex; j < endIndex; j++)
{ {
s += values[j].Magnitude; if (i == aat.ColumnIndices[j])
{
norm += aat.Values[j].Magnitude;
}
} }
norm = Math.Max(norm, s);
} }
return Math.Sqrt(norm);
return norm;
} }
/// <summary> /// <summary>

14
src/Numerics/LinearAlgebra/Double/DenseMatrix.cs

@ -425,13 +425,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double
return Control.LinearAlgebraProvider.MatrixNorm(Norm.OneNorm, _rowCount, _columnCount, _values); return Control.LinearAlgebraProvider.MatrixNorm(Norm.OneNorm, _rowCount, _columnCount, _values);
} }
/// <summary>Calculates the Frobenius norm of this matrix.</summary>
/// <returns>The Frobenius norm of this matrix.</returns>
public override double FrobeniusNorm()
{
return Control.LinearAlgebraProvider.MatrixNorm(Norm.FrobeniusNorm, _rowCount, _columnCount, _values);
}
/// <summary>Calculates the infinity norm of this matrix.</summary> /// <summary>Calculates the infinity norm of this matrix.</summary>
/// <returns>The infinity norm of this matrix.</returns> /// <returns>The infinity norm of this matrix.</returns>
public override double InfinityNorm() public override double InfinityNorm()
@ -439,6 +432,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double
return Control.LinearAlgebraProvider.MatrixNorm(Norm.InfinityNorm, _rowCount, _columnCount, _values); return Control.LinearAlgebraProvider.MatrixNorm(Norm.InfinityNorm, _rowCount, _columnCount, _values);
} }
/// <summary>Calculates the Frobenius norm of this matrix.</summary>
/// <returns>The Frobenius norm of this matrix.</returns>
public override double FrobeniusNorm()
{
return Control.LinearAlgebraProvider.MatrixNorm(Norm.FrobeniusNorm, _rowCount, _columnCount, _values);
}
/// <summary> /// <summary>
/// Add a scalar to each element of the matrix and stores the result in the result vector. /// Add a scalar to each element of the matrix and stores the result in the result vector.
/// </summary> /// </summary>

15
src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs

@ -655,14 +655,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double
return _data.Aggregate(double.NegativeInfinity, (current, t) => Math.Max(current, Math.Abs(t))); return _data.Aggregate(double.NegativeInfinity, (current, t) => Math.Max(current, Math.Abs(t)));
} }
/// <summary>Calculates the Frobenius norm of this matrix.</summary>
/// <returns>The Frobenius norm of this matrix.</returns>
public override double FrobeniusNorm()
{
var norm = _data.Sum(t => t * t);
return Math.Sqrt(norm);
}
/// <summary>Calculates the infinity norm of this matrix.</summary> /// <summary>Calculates the infinity norm of this matrix.</summary>
/// <returns>The infinity norm of this matrix.</returns> /// <returns>The infinity norm of this matrix.</returns>
public override double InfinityNorm() public override double InfinityNorm()
@ -670,6 +662,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double
return L1Norm(); return L1Norm();
} }
/// <summary>Calculates the Frobenius norm of this matrix.</summary>
/// <returns>The Frobenius norm of this matrix.</returns>
public override double FrobeniusNorm()
{
return Math.Sqrt(_data.Sum(t => t * t));
}
/// <summary>Calculates the condition number of this matrix.</summary> /// <summary>Calculates the condition number of this matrix.</summary>
/// <returns>The condition number of the matrix.</returns> /// <returns>The condition number of the matrix.</returns>
public override double ConditionNumber() public override double ConditionNumber()

2
src/Numerics/LinearAlgebra/Double/Factorization/Svd.cs

@ -72,7 +72,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
/// Gets the two norm of the <see cref="Matrix{T}"/>. /// Gets the two norm of the <see cref="Matrix{T}"/>.
/// </summary> /// </summary>
/// <returns>The 2-norm of the <see cref="Matrix{T}"/>.</returns> /// <returns>The 2-norm of the <see cref="Matrix{T}"/>.</returns>
public override double Norm2 public override double L2Norm
{ {
get get
{ {

56
src/Numerics/LinearAlgebra/Double/Matrix.cs

@ -54,28 +54,34 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <returns>The L1 norm of the matrix.</returns> /// <returns>The L1 norm of the matrix.</returns>
public override double L1Norm() public override double L1Norm()
{ {
var norm = 0.0; var norm = 0d;
for (var j = 0; j < ColumnCount; j++) for (var j = 0; j < ColumnCount; j++)
{ {
var s = 0.0; var s = 0d;
for (var i = 0; i < RowCount; i++) for (var i = 0; i < RowCount; i++)
{ {
s += Math.Abs(At(i, j)); s += Math.Abs(At(i, j));
} }
norm = Math.Max(norm, s); norm = Math.Max(norm, s);
} }
return norm; return norm;
} }
/// <summary> /// <summary>Calculates the infinity norm of this matrix.</summary>
/// Returns the conjugate transpose of this matrix. /// <returns>The infinity norm of this matrix.</returns>
/// </summary> public override double InfinityNorm()
/// <returns>The conjugate transpose of this matrix.</returns>
public override sealed Matrix<double> ConjugateTranspose()
{ {
return Transpose(); var norm = 0d;
for (var i = 0; i < RowCount; i++)
{
var s = 0d;
for (var j = 0; j < ColumnCount; j++)
{
s += Math.Abs(At(i, j));
}
norm = Math.Max(norm, s);
}
return norm;
} }
/// <summary>Calculates the Frobenius norm of this matrix.</summary> /// <summary>Calculates the Frobenius norm of this matrix.</summary>
@ -84,35 +90,21 @@ namespace MathNet.Numerics.LinearAlgebra.Double
{ {
var transpose = Transpose(); var transpose = Transpose();
var aat = this * transpose; var aat = this * transpose;
var norm = 0d;
var norm = 0.0;
for (var i = 0; i < RowCount; i++) for (var i = 0; i < RowCount; i++)
{ {
norm += aat.At(i, i); norm += aat.At(i, i);
} }
return Math.Sqrt(norm);
norm = Math.Sqrt(norm);
return norm;
} }
/// <summary>Calculates the infinity norm of this matrix.</summary> /// <summary>
/// <returns>The infinity norm of this matrix.</returns> /// Returns the conjugate transpose of this matrix.
public override double InfinityNorm() /// </summary>
/// <returns>The conjugate transpose of this matrix.</returns>
public override sealed Matrix<double> ConjugateTranspose()
{ {
var norm = 0.0; return Transpose();
for (var i = 0; i < RowCount; i++)
{
var s = 0.0;
for (var j = 0; j < ColumnCount; j++)
{
s += Math.Abs(At(i, j));
}
norm = Math.Max(norm, s);
}
return norm;
} }
/// <summary> /// <summary>

53
src/Numerics/LinearAlgebra/Double/SparseMatrix.cs

@ -668,17 +668,17 @@ namespace MathNet.Numerics.LinearAlgebra.Double
return new SparseMatrix(ret); return new SparseMatrix(ret);
} }
/// <summary>Calculates the Frobenius norm of this matrix.</summary> /// <summary>Calculates the infinity norm of this matrix.</summary>
/// <returns>The Frobenius norm of this matrix.</returns> /// <returns>The infinity norm of this matrix.</returns>
public override double FrobeniusNorm() public override double InfinityNorm()
{ {
var aat = (SparseCompressedRowMatrixStorage<double>) (this*Transpose()).Storage; var rowPointers = _storage.RowPointers;
var values = _storage.Values;
var norm = 0d; var norm = 0d;
for (var i = 0; i < RowCount; i++)
for (var i = 0; i < aat.RowCount; i++)
{ {
var startIndex = aat.RowPointers[i]; var startIndex = rowPointers[i];
var endIndex = aat.RowPointers[i + 1]; var endIndex = rowPointers[i + 1];
if (startIndex == endIndex) if (startIndex == endIndex)
{ {
@ -686,30 +686,26 @@ namespace MathNet.Numerics.LinearAlgebra.Double
continue; continue;
} }
var s = 0d;
for (var j = startIndex; j < endIndex; j++) for (var j = startIndex; j < endIndex; j++)
{ {
if (i == aat.ColumnIndices[j]) s += Math.Abs(values[j]);
{
norm += Math.Abs(aat.Values[j]);
}
} }
norm = Math.Max(norm, s);
} }
return norm;
return Math.Sqrt(norm);
} }
/// <summary>Calculates the infinity norm of this matrix.</summary> /// <summary>Calculates the Frobenius norm of this matrix.</summary>
/// <returns>The infinity norm of this matrix.</returns> /// <returns>The Frobenius norm of this matrix.</returns>
public override double InfinityNorm() public override double FrobeniusNorm()
{ {
var rowPointers = _storage.RowPointers; var aat = (SparseCompressedRowMatrixStorage<double>) (this*Transpose()).Storage;
var values = _storage.Values;
var norm = 0d; var norm = 0d;
for (var i = 0; i < RowCount; i++) for (var i = 0; i < aat.RowCount; i++)
{ {
var startIndex = rowPointers[i]; var startIndex = aat.RowPointers[i];
var endIndex = rowPointers[i + 1]; var endIndex = aat.RowPointers[i + 1];
if (startIndex == endIndex) if (startIndex == endIndex)
{ {
@ -717,16 +713,15 @@ namespace MathNet.Numerics.LinearAlgebra.Double
continue; continue;
} }
var s = 0d;
for (var j = startIndex; j < endIndex; j++) for (var j = startIndex; j < endIndex; j++)
{ {
s += Math.Abs(values[j]); if (i == aat.ColumnIndices[j])
{
norm += Math.Abs(aat.Values[j]);
}
} }
norm = Math.Max(norm, s);
} }
return Math.Sqrt(norm);
return norm;
} }
/// <summary> /// <summary>

2
src/Numerics/LinearAlgebra/Factorization/Svd.cs

@ -120,7 +120,7 @@ namespace MathNet.Numerics.LinearAlgebra.Factorization
/// Gets the two norm of the <see cref="Matrix{T}"/>. /// Gets the two norm of the <see cref="Matrix{T}"/>.
/// </summary> /// </summary>
/// <returns>The 2-norm of the <see cref="Matrix{T}"/>.</returns> /// <returns>The 2-norm of the <see cref="Matrix{T}"/>.</returns>
public abstract T Norm2 { get; } public abstract double L2Norm { get; }
/// <summary> /// <summary>
/// Gets the condition number <b>max(S) / min(S)</b> /// Gets the condition number <b>max(S) / min(S)</b>

21
src/Numerics/LinearAlgebra/Matrix.Arithmetic.cs

@ -1148,6 +1148,27 @@ namespace MathNet.Numerics.LinearAlgebra
return ret; return ret;
} }
/// <summary>Calculates the L1 norm.</summary>
/// <returns>The L1 norm of the matrix.</returns>
public abstract double L1Norm();
/// <summary>Calculates the L2 norm.</summary>
/// <returns>The L2 norm of the matrix.</returns>
/// <remarks>For sparse matrices, the L2 norm is computed using a dense implementation of singular value decomposition.
/// In a later release, it will be replaced with a sparse implementation.</remarks>
public virtual double L2Norm()
{
return Svd(false).L2Norm;
}
/// <summary>Calculates the infinity norm of this matrix.</summary>
/// <returns>The infinity norm of this matrix.</returns>
public abstract double InfinityNorm();
/// <summary>Calculates the Frobenius norm of this matrix.</summary>
/// <returns>The Frobenius norm of this matrix.</returns>
public abstract double FrobeniusNorm();
#region Exceptions - possibly move elsewhere? #region Exceptions - possibly move elsewhere?
internal static Exception DimensionsDontMatch<TException>(Matrix<T> left, Matrix<T> right, Matrix<T> result, string paramName = null) internal static Exception DimensionsDontMatch<TException>(Matrix<T> left, Matrix<T> right, Matrix<T> result, string paramName = null)

21
src/Numerics/LinearAlgebra/Matrix.cs

@ -1251,27 +1251,6 @@ namespace MathNet.Numerics.LinearAlgebra
lower.Storage.CopySubMatrixToUnchecked(result.Storage, 0, RowCount, lower.RowCount, 0, ColumnCount, lower.ColumnCount); lower.Storage.CopySubMatrixToUnchecked(result.Storage, 0, RowCount, lower.RowCount, 0, ColumnCount, lower.ColumnCount);
} }
/// <summary>Calculates the L1 norm.</summary>
/// <returns>The L1 norm of the matrix.</returns>
public abstract T L1Norm();
/// <summary>Calculates the L2 norm.</summary>
/// <returns>The L2 norm of the matrix.</returns>
/// <remarks>For sparse matrices, the L2 norm is computed using a dense implementation of singular value decomposition.
/// In a later release, it will be replaced with a sparse implementation.</remarks>
public virtual T L2Norm()
{
return Svd(false).Norm2;
}
/// <summary>Calculates the Frobenius norm of this matrix.</summary>
/// <returns>The Frobenius norm of this matrix.</returns>
public abstract T FrobeniusNorm();
/// <summary>Calculates the infinity norm of this matrix.</summary>
/// <returns>The infinity norm of this matrix.</returns>
public abstract T InfinityNorm();
/// <summary> /// <summary>
/// Gets a value indicating whether this matrix is symmetric. /// Gets a value indicating whether this matrix is symmetric.
/// </summary> /// </summary>

18
src/Numerics/LinearAlgebra/Single/DenseMatrix.cs

@ -420,25 +420,25 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <summary>Calculates the L1 norm.</summary> /// <summary>Calculates the L1 norm.</summary>
/// <returns>The L1 norm of the matrix.</returns> /// <returns>The L1 norm of the matrix.</returns>
public override float L1Norm() public override double L1Norm()
{ {
return Control.LinearAlgebraProvider.MatrixNorm(Norm.OneNorm, _rowCount, _columnCount, _values); return Control.LinearAlgebraProvider.MatrixNorm(Norm.OneNorm, _rowCount, _columnCount, _values);
} }
/// <summary>Calculates the Frobenius norm of this matrix.</summary>
/// <returns>The Frobenius norm of this matrix.</returns>
public override float FrobeniusNorm()
{
return Control.LinearAlgebraProvider.MatrixNorm(Norm.FrobeniusNorm, _rowCount, _columnCount, _values);
}
/// <summary>Calculates the infinity norm of this matrix.</summary> /// <summary>Calculates the infinity norm of this matrix.</summary>
/// <returns>The infinity norm of this matrix.</returns> /// <returns>The infinity norm of this matrix.</returns>
public override float InfinityNorm() public override double InfinityNorm()
{ {
return Control.LinearAlgebraProvider.MatrixNorm(Norm.InfinityNorm, _rowCount, _columnCount, _values); return Control.LinearAlgebraProvider.MatrixNorm(Norm.InfinityNorm, _rowCount, _columnCount, _values);
} }
/// <summary>Calculates the Frobenius norm of this matrix.</summary>
/// <returns>The Frobenius norm of this matrix.</returns>
public override double FrobeniusNorm()
{
return Control.LinearAlgebraProvider.MatrixNorm(Norm.FrobeniusNorm, _rowCount, _columnCount, _values);
}
/// <summary> /// <summary>
/// Add a scalar to each element of the matrix and stores the result in the result vector. /// Add a scalar to each element of the matrix and stores the result in the result vector.
/// </summary> /// </summary>

21
src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs

@ -643,33 +643,32 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <summary>Calculates the L1 norm.</summary> /// <summary>Calculates the L1 norm.</summary>
/// <returns>The L1 norm of the matrix.</returns> /// <returns>The L1 norm of the matrix.</returns>
public override float L1Norm() public override double L1Norm()
{ {
return _data.Aggregate(float.NegativeInfinity, (current, t) => Math.Max(current, Math.Abs(t))); return _data.Aggregate(float.NegativeInfinity, (current, t) => Math.Max(current, Math.Abs(t)));
} }
/// <summary>Calculates the L2 norm.</summary> /// <summary>Calculates the L2 norm.</summary>
/// <returns>The L2 norm of the matrix.</returns> /// <returns>The L2 norm of the matrix.</returns>
public override float L2Norm() public override double L2Norm()
{ {
return _data.Aggregate(float.NegativeInfinity, (current, t) => Math.Max(current, Math.Abs(t))); return _data.Aggregate(float.NegativeInfinity, (current, t) => Math.Max(current, Math.Abs(t)));
} }
/// <summary>Calculates the Frobenius norm of this matrix.</summary>
/// <returns>The Frobenius norm of this matrix.</returns>
public override float FrobeniusNorm()
{
var norm = _data.Sum(t => t * t);
return Convert.ToSingle(Math.Sqrt(norm));
}
/// <summary>Calculates the infinity norm of this matrix.</summary> /// <summary>Calculates the infinity norm of this matrix.</summary>
/// <returns>The infinity norm of this matrix.</returns> /// <returns>The infinity norm of this matrix.</returns>
public override float InfinityNorm() public override double InfinityNorm()
{ {
return L1Norm(); return L1Norm();
} }
/// <summary>Calculates the Frobenius norm of this matrix.</summary>
/// <returns>The Frobenius norm of this matrix.</returns>
public override double FrobeniusNorm()
{
return Math.Sqrt(_data.Sum(t => t * t));
}
/// <summary>Calculates the condition number of this matrix.</summary> /// <summary>Calculates the condition number of this matrix.</summary>
/// <returns>The condition number of the matrix.</returns> /// <returns>The condition number of the matrix.</returns>
public override float ConditionNumber() public override float ConditionNumber()

2
src/Numerics/LinearAlgebra/Single/Factorization/Svd.cs

@ -72,7 +72,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
/// Gets the two norm of the <see cref="Matrix{T}"/>. /// Gets the two norm of the <see cref="Matrix{T}"/>.
/// </summary> /// </summary>
/// <returns>The 2-norm of the <see cref="Matrix{T}"/>.</returns> /// <returns>The 2-norm of the <see cref="Matrix{T}"/>.</returns>
public override float Norm2 public override double L2Norm
{ {
get get
{ {

60
src/Numerics/LinearAlgebra/Single/Matrix.cs

@ -52,67 +52,59 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <summary>Calculates the L1 norm.</summary> /// <summary>Calculates the L1 norm.</summary>
/// <returns>The L1 norm of the matrix.</returns> /// <returns>The L1 norm of the matrix.</returns>
public override float L1Norm() public override double L1Norm()
{ {
var norm = 0.0f; var norm = 0d;
for (var j = 0; j < ColumnCount; j++) for (var j = 0; j < ColumnCount; j++)
{ {
var s = 0.0f; var s = 0d;
for (var i = 0; i < RowCount; i++) for (var i = 0; i < RowCount; i++)
{ {
s += Math.Abs(At(i, j)); s += Math.Abs(At(i, j));
} }
norm = Math.Max(norm, s); norm = Math.Max(norm, s);
} }
return norm; return norm;
} }
/// <summary> /// <summary>Calculates the infinity norm of this matrix.</summary>
/// Returns the conjugate transpose of this matrix. /// <returns>The infinity norm of this matrix.</returns>
/// </summary> public override double InfinityNorm()
/// <returns>The conjugate transpose of this matrix.</returns>
public override sealed Matrix<float> ConjugateTranspose()
{ {
return Transpose(); var norm = 0d;
for (var i = 0; i < RowCount; i++)
{
var s = 0d;
for (var j = 0; j < ColumnCount; j++)
{
s += Math.Abs(At(i, j));
}
norm = Math.Max(norm, s);
}
return norm;
} }
/// <summary>Calculates the Frobenius norm of this matrix.</summary> /// <summary>Calculates the Frobenius norm of this matrix.</summary>
/// <returns>The Frobenius norm of this matrix.</returns> /// <returns>The Frobenius norm of this matrix.</returns>
public override float FrobeniusNorm() public override double FrobeniusNorm()
{ {
var transpose = Transpose(); var transpose = Transpose();
var aat = this * transpose; var aat = this * transpose;
var norm = 0d;
var norm = 0.0f;
for (var i = 0; i < RowCount; i++) for (var i = 0; i < RowCount; i++)
{ {
norm += aat.At(i, i); norm += aat.At(i, i);
} }
return Math.Sqrt(norm);
norm = Convert.ToSingle(Math.Sqrt(norm));
return norm;
} }
/// <summary>Calculates the infinity norm of this matrix.</summary> /// <summary>
/// <returns>The infinity norm of this matrix.</returns> /// Returns the conjugate transpose of this matrix.
public override float InfinityNorm() /// </summary>
/// <returns>The conjugate transpose of this matrix.</returns>
public override sealed Matrix<float> ConjugateTranspose()
{ {
var norm = 0.0f; return Transpose();
for (var i = 0; i < RowCount; i++)
{
var s = 0.0f;
for (var j = 0; j < ColumnCount; j++)
{
s += Math.Abs(At(i, j));
}
norm = Math.Max(norm, s);
}
return norm;
} }
/// <summary> /// <summary>

57
src/Numerics/LinearAlgebra/Single/SparseMatrix.cs

@ -670,18 +670,18 @@ namespace MathNet.Numerics.LinearAlgebra.Single
return new SparseMatrix(ret); return new SparseMatrix(ret);
} }
/// <summary>Calculates the Frobenius norm of this matrix.</summary> /// <summary>Calculates the infinity norm of this matrix.</summary>
/// <returns>The Frobenius norm of this matrix.</returns> /// <returns>The infinity norm of this matrix.</returns>
public override float FrobeniusNorm() public override double InfinityNorm()
{ {
var aat = (SparseCompressedRowMatrixStorage<float>) (this*Transpose()).Storage; var rowPointers = _storage.RowPointers;
var norm = 0f; var values = _storage.Values;
var norm = 0d;
for (var i = 0; i < aat.RowCount; i++) for (var i = 0; i < RowCount; i++)
{ {
// Get the begin / end index for the current row // Get the begin / end index for the current row
var startIndex = aat.RowPointers[i]; var startIndex = rowPointers[i];
var endIndex = aat.RowPointers[i + 1]; var endIndex = rowPointers[i + 1];
// Get the values for the current row // Get the values for the current row
if (startIndex == endIndex) if (startIndex == endIndex)
@ -690,31 +690,27 @@ namespace MathNet.Numerics.LinearAlgebra.Single
continue; continue;
} }
var s = 0d;
for (var j = startIndex; j < endIndex; j++) for (var j = startIndex; j < endIndex; j++)
{ {
if (i == aat.ColumnIndices[j]) s += Math.Abs(values[j]);
{
norm += Math.Abs(aat.Values[j]);
}
} }
norm = Math.Max(norm, s);
} }
return norm;
return Convert.ToSingle(Math.Sqrt(norm));
} }
/// <summary>Calculates the infinity norm of this matrix.</summary> /// <summary>Calculates the Frobenius norm of this matrix.</summary>
/// <returns>The infinity norm of this matrix.</returns> /// <returns>The Frobenius norm of this matrix.</returns>
public override float InfinityNorm() public override double FrobeniusNorm()
{ {
var rowPointers = _storage.RowPointers; var aat = (SparseCompressedRowMatrixStorage<float>) (this*Transpose()).Storage;
var values = _storage.Values; var norm = 0d;
for (var i = 0; i < aat.RowCount; i++)
var norm = 0f;
for (var i = 0; i < RowCount; i++)
{ {
// Get the begin / end index for the current row // Get the begin / end index for the current row
var startIndex = rowPointers[i]; var startIndex = aat.RowPointers[i];
var endIndex = rowPointers[i + 1]; var endIndex = aat.RowPointers[i + 1];
// Get the values for the current row // Get the values for the current row
if (startIndex == endIndex) if (startIndex == endIndex)
@ -723,16 +719,15 @@ namespace MathNet.Numerics.LinearAlgebra.Single
continue; continue;
} }
var s = 0f;
for (var j = startIndex; j < endIndex; j++) for (var j = startIndex; j < endIndex; j++)
{ {
s += Math.Abs(values[j]); if (i == aat.ColumnIndices[j])
{
norm += Math.Abs(aat.Values[j]);
}
} }
norm = Math.Max(norm, s);
} }
return Math.Sqrt(norm);
return norm;
} }
/// <summary> /// <summary>

138
src/Numerics/Precision.Equality.cs

@ -30,10 +30,11 @@
using System; using System;
using System.Collections.Generic; using System.Collections.Generic;
using MathNet.Numerics.LinearAlgebra;
#if PORTABLE #if PORTABLE
using System.Runtime.InteropServices; using System.Runtime.InteropServices;
#endif #endif
using MathNet.Numerics.LinearAlgebra;
namespace MathNet.Numerics namespace MathNet.Numerics
{ {
@ -42,6 +43,8 @@ namespace MathNet.Numerics
using Complex = System.Numerics.Complex; using Complex = System.Numerics.Complex;
#endif #endif
// TODO PERF: Cache/Precompute 10^x terms
public static partial class Precision public static partial class Precision
{ {
/// <summary> /// <summary>
@ -116,7 +119,7 @@ namespace MathNet.Numerics
// If one is almost zero, fall back to absolute equality // If one is almost zero, fall back to absolute equality
if (Math.Abs(a) < DoublePrecision || Math.Abs(b) < DoublePrecision) if (Math.Abs(a) < DoublePrecision || Math.Abs(b) < DoublePrecision)
{ {
return AlmostEqualNorm(a, b, diff, maximumError); return Math.Abs(diff) < maximumError;
} }
if ((a == 0 && Math.Abs(b) < maximumError) || (b == 0 && Math.Abs(a) < maximumError)) if ((a == 0 && Math.Abs(b) < maximumError) || (b == 0 && Math.Abs(a) < maximumError))
@ -361,8 +364,7 @@ namespace MathNet.Numerics
// 10^(-numberOfDecimalPlaces). We divide by two so that we have half the range // 10^(-numberOfDecimalPlaces). We divide by two so that we have half the range
// on each side of the numbers, e.g. if decimalPlaces == 2, // on each side of the numbers, e.g. if decimalPlaces == 2,
// then 0.01 will equal between 0.005 and 0.015, but not 0.02 and not 0.00 // then 0.01 will equal between 0.005 and 0.015, but not 0.02 and not 0.00
double decimalPlaceMagnitude = Math.Pow(10, -decimalPlaces) / 2d; return Math.Abs(diff) < Math.Pow(10, -decimalPlaces) / 2d;
return Math.Abs(diff) < decimalPlaceMagnitude;
} }
/// <summary> /// <summary>
@ -437,8 +439,7 @@ namespace MathNet.Numerics
// 10^(-numberOfDecimalPlaces). We divide by two so that we have half the range // 10^(-numberOfDecimalPlaces). We divide by two so that we have half the range
// on each side of the numbers, e.g. if decimalPlaces == 2, // on each side of the numbers, e.g. if decimalPlaces == 2,
// then 0.01 will equal between 0.005 and 0.015, but not 0.02 and not 0.00 // then 0.01 will equal between 0.005 and 0.015, but not 0.02 and not 0.00
double decimalPlaceMagnitude = Math.Pow(10, -decimalPlaces) / 2d; return Math.Abs(diff) < Math.Pow(10, -decimalPlaces) / 2d;
return Math.Abs(diff) < decimalPlaceMagnitude;
} }
// If the magnitudes of the two numbers are equal to within one magnitude the numbers could potentially be equal // If the magnitudes of the two numbers are equal to within one magnitude the numbers could potentially be equal
@ -454,8 +455,7 @@ namespace MathNet.Numerics
// 10^(-numberOfDecimalPlaces). We divide by two so that we have half the range // 10^(-numberOfDecimalPlaces). We divide by two so that we have half the range
// on each side of the numbers, e.g. if decimalPlaces == 2, // on each side of the numbers, e.g. if decimalPlaces == 2,
// then 0.01 will equal between 0.00995 and 0.01005, but not 0.0015 and not 0.0095 // then 0.01 will equal between 0.00995 and 0.01005, but not 0.0015 and not 0.0095
double decimalPlaceMagnitude1 = Math.Pow(10, magnitudeOfMax - decimalPlaces) / 2d; return Math.Abs(diff) < Math.Pow(10, magnitudeOfMax - decimalPlaces) / 2d;
return Math.Abs(diff) < decimalPlaceMagnitude1;
} }
/// <summary> /// <summary>
@ -775,6 +775,7 @@ namespace MathNet.Numerics
/// </summary> /// </summary>
/// <param name="a">The first value list.</param> /// <param name="a">The first value list.</param>
/// <param name="b">The second value list.</param> /// <param name="b">The second value list.</param>
/// <param name="decimalPlaces">The number of decimal places.</param>
public static bool ListAlmostEqual(this IList<double> a, IList<double> b, int decimalPlaces) public static bool ListAlmostEqual(this IList<double> a, IList<double> b, int decimalPlaces)
{ {
return ListForAll(a, b, AlmostEqual, decimalPlaces); return ListForAll(a, b, AlmostEqual, decimalPlaces);
@ -786,6 +787,7 @@ namespace MathNet.Numerics
/// </summary> /// </summary>
/// <param name="a">The first value list.</param> /// <param name="a">The first value list.</param>
/// <param name="b">The second value list.</param> /// <param name="b">The second value list.</param>
/// <param name="decimalPlaces">The number of decimal places.</param>
public static bool ListAlmostEqual(this IList<float> a, IList<float> b, int decimalPlaces) public static bool ListAlmostEqual(this IList<float> a, IList<float> b, int decimalPlaces)
{ {
return ListForAll(a, b, AlmostEqual, decimalPlaces); return ListForAll(a, b, AlmostEqual, decimalPlaces);
@ -797,6 +799,7 @@ namespace MathNet.Numerics
/// </summary> /// </summary>
/// <param name="a">The first value list.</param> /// <param name="a">The first value list.</param>
/// <param name="b">The second value list.</param> /// <param name="b">The second value list.</param>
/// <param name="decimalPlaces">The number of decimal places.</param>
public static bool ListAlmostEqual(this IList<Complex> a, IList<Complex> b, int decimalPlaces) public static bool ListAlmostEqual(this IList<Complex> a, IList<Complex> b, int decimalPlaces)
{ {
return ListForAll(a, b, AlmostEqual, decimalPlaces); return ListForAll(a, b, AlmostEqual, decimalPlaces);
@ -808,6 +811,7 @@ namespace MathNet.Numerics
/// </summary> /// </summary>
/// <param name="a">The first value list.</param> /// <param name="a">The first value list.</param>
/// <param name="b">The second value list.</param> /// <param name="b">The second value list.</param>
/// <param name="decimalPlaces">The number of decimal places.</param>
public static bool ListAlmostEqual(this IList<Complex32> a, IList<Complex32> b, int decimalPlaces) public static bool ListAlmostEqual(this IList<Complex32> a, IList<Complex32> b, int decimalPlaces)
{ {
return ListForAll(a, b, AlmostEqual, decimalPlaces); return ListForAll(a, b, AlmostEqual, decimalPlaces);
@ -819,6 +823,7 @@ namespace MathNet.Numerics
/// </summary> /// </summary>
/// <param name="a">The first value list.</param> /// <param name="a">The first value list.</param>
/// <param name="b">The second value list.</param> /// <param name="b">The second value list.</param>
/// <param name="decimalPlaces">The number of decimal places.</param>
public static bool ListAlmostEqualRelative(this IList<double> a, IList<double> b, int decimalPlaces) public static bool ListAlmostEqualRelative(this IList<double> a, IList<double> b, int decimalPlaces)
{ {
return ListForAll(a, b, AlmostEqualRelative, decimalPlaces); return ListForAll(a, b, AlmostEqualRelative, decimalPlaces);
@ -830,6 +835,7 @@ namespace MathNet.Numerics
/// </summary> /// </summary>
/// <param name="a">The first value list.</param> /// <param name="a">The first value list.</param>
/// <param name="b">The second value list.</param> /// <param name="b">The second value list.</param>
/// <param name="decimalPlaces">The number of decimal places.</param>
public static bool ListAlmostEqualRelative(this IList<float> a, IList<float> b, int decimalPlaces) public static bool ListAlmostEqualRelative(this IList<float> a, IList<float> b, int decimalPlaces)
{ {
return ListForAll(a, b, AlmostEqualRelative, decimalPlaces); return ListForAll(a, b, AlmostEqualRelative, decimalPlaces);
@ -841,6 +847,7 @@ namespace MathNet.Numerics
/// </summary> /// </summary>
/// <param name="a">The first value list.</param> /// <param name="a">The first value list.</param>
/// <param name="b">The second value list.</param> /// <param name="b">The second value list.</param>
/// <param name="decimalPlaces">The number of decimal places.</param>
public static bool ListAlmostEqualRelative(this IList<Complex> a, IList<Complex> b, int decimalPlaces) public static bool ListAlmostEqualRelative(this IList<Complex> a, IList<Complex> b, int decimalPlaces)
{ {
return ListForAll(a, b, AlmostEqualRelative, decimalPlaces); return ListForAll(a, b, AlmostEqualRelative, decimalPlaces);
@ -852,6 +859,7 @@ namespace MathNet.Numerics
/// </summary> /// </summary>
/// <param name="a">The first value list.</param> /// <param name="a">The first value list.</param>
/// <param name="b">The second value list.</param> /// <param name="b">The second value list.</param>
/// <param name="decimalPlaces">The number of decimal places.</param>
public static bool ListAlmostEqualRelative(this IList<Complex32> a, IList<Complex32> b, int decimalPlaces) public static bool ListAlmostEqualRelative(this IList<Complex32> a, IList<Complex32> b, int decimalPlaces)
{ {
return ListForAll(a, b, AlmostEqualRelative, decimalPlaces); return ListForAll(a, b, AlmostEqualRelative, decimalPlaces);
@ -943,8 +951,7 @@ namespace MathNet.Numerics
} }
/// <summary> /// <summary>
/// Compares two complex and determines if they are equal within /// Compares two vectors and determines if they are equal within the specified maximum error.
/// the specified maximum error.
/// </summary> /// </summary>
/// <param name="a">The first value.</param> /// <param name="a">The first value.</param>
/// <param name="b">The second value.</param> /// <param name="b">The second value.</param>
@ -956,8 +963,7 @@ namespace MathNet.Numerics
} }
/// <summary> /// <summary>
/// Compares two doubles and determines if they are equal within /// Compares two vectors and determines if they are equal within the specified maximum error.
/// the specified maximum error.
/// </summary> /// </summary>
/// <param name="a">The first value.</param> /// <param name="a">The first value.</param>
/// <param name="b">The second value.</param> /// <param name="b">The second value.</param>
@ -969,8 +975,8 @@ namespace MathNet.Numerics
} }
/// <summary> /// <summary>
/// Compares two doubles and determines if they are equal to within the specified number of decimal places or not, using the /// Compares two vectors and determines if they are equal to within the specified number
/// number of decimal places as an absolute measure. /// of decimal places or not, using the number of decimal places as an absolute measure.
/// </summary> /// </summary>
/// <param name="a">The first value.</param> /// <param name="a">The first value.</param>
/// <param name="b">The second value.</param> /// <param name="b">The second value.</param>
@ -982,8 +988,8 @@ namespace MathNet.Numerics
} }
/// <summary> /// <summary>
/// Compares two doubles and determines if they are equal to within the specified number of decimal places or not. If the numbers /// Compares two vectors and determines if they are equal to within the specified number of decimal places or not.
/// are very close to zero an absolute difference is compared, otherwise the relative difference is compared. /// If the numbers are very close to zero an absolute difference is compared, otherwise the relative difference is compared.
/// </summary> /// </summary>
/// <param name="a">The first value.</param> /// <param name="a">The first value.</param>
/// <param name="b">The second value.</param> /// <param name="b">The second value.</param>
@ -994,56 +1000,54 @@ namespace MathNet.Numerics
return AlmostEqualNormRelative(a.L2Norm(), b.L2Norm(), (a - b).L2Norm(), decimalPlaces); return AlmostEqualNormRelative(a.L2Norm(), b.L2Norm(), (a - b).L2Norm(), decimalPlaces);
} }
///// <summary> /// <summary>
///// Compares two complex and determines if they are equal within /// Compares two matrices and determines if they are equal within the specified maximum error.
///// the specified maximum error. /// </summary>
///// </summary> /// <param name="a">The first value.</param>
///// <param name="a">The first value.</param> /// <param name="b">The second value.</param>
///// <param name="b">The second value.</param> /// <param name="maximumAbsoluteError">The accuracy required for being almost equal.</param>
///// <param name="maximumAbsoluteError">The accuracy required for being almost equal.</param> public static bool AlmostEqual<T>(this Matrix<T> a, Matrix<T> b, double maximumAbsoluteError)
//public static bool AlmostEqual<T>(this Matrix<T> a, Matrix<T> b, double maximumAbsoluteError) where T : struct, IEquatable<T>, IFormattable
// where T : struct, IEquatable<T>, IFormattable {
//{ return AlmostEqualNorm(a.L2Norm(), b.L2Norm(), (a - b).L2Norm(), maximumAbsoluteError);
// return AlmostEqualNorm(a.L2Norm(), b.L2Norm(), (a - b).L2Norm(), maximumAbsoluteError); }
//}
/// <summary>
///// <summary> /// Compares two matrices and determines if they are equal within the specified maximum error.
///// Compares two doubles and determines if they are equal within /// </summary>
///// the specified maximum error. /// <param name="a">The first value.</param>
///// </summary> /// <param name="b">The second value.</param>
///// <param name="a">The first value.</param> /// <param name="maximumError">The accuracy required for being almost equal.</param>
///// <param name="b">The second value.</param> public static bool AlmostEqualRelative<T>(this Matrix<T> a, Matrix<T> b, double maximumError)
///// <param name="maximumError">The accuracy required for being almost equal.</param> where T : struct, IEquatable<T>, IFormattable
//public static bool AlmostEqualRelative<T>(this Matrix<T> a, Matrix<T> b, double maximumError) {
// where T : struct, IEquatable<T>, IFormattable return AlmostEqualNormRelative(a.L2Norm(), b.L2Norm(), (a - b).L2Norm(), maximumError);
//{ }
// return AlmostEqualNormRelative(a.L2Norm(), b.L2Norm(), (a - b).L2Norm(), maximumError);
//} /// <summary>
/// Compares two matrices and determines if they are equal to within the specified number
///// <summary> /// of decimal places or not, using the number of decimal places as an absolute measure.
///// Compares two doubles and determines if they are equal to within the specified number of decimal places or not, using the /// </summary>
///// number of decimal places as an absolute measure. /// <param name="a">The first value.</param>
///// </summary> /// <param name="b">The second value.</param>
///// <param name="a">The first value.</param> /// <param name="decimalPlaces">The number of decimal places.</param>
///// <param name="b">The second value.</param> public static bool AlmostEqual<T>(this Matrix<T> a, Matrix<T> b, int decimalPlaces)
///// <param name="decimalPlaces">The number of decimal places.</param> where T : struct, IEquatable<T>, IFormattable
//public static bool AlmostEqual<T>(this Matrix<T> a, Matrix<T> b, int decimalPlaces) {
// where T : struct, IEquatable<T>, IFormattable return AlmostEqualNorm(a.L2Norm(), b.L2Norm(), (a - b).L2Norm(), decimalPlaces);
//{ }
// return AlmostEqualNorm(a.L2Norm(), b.L2Norm(), (a - b).L2Norm(), decimalPlaces);
//} /// <summary>
/// Compares two matrices and determines if they are equal to within the specified number of decimal places or not.
///// <summary> /// If the numbers are very close to zero an absolute difference is compared, otherwise the relative difference is compared.
///// Compares two doubles and determines if they are equal to within the specified number of decimal places or not. If the numbers /// </summary>
///// are very close to zero an absolute difference is compared, otherwise the relative difference is compared. /// <param name="a">The first value.</param>
///// </summary> /// <param name="b">The second value.</param>
///// <param name="a">The first value.</param> /// <param name="decimalPlaces">The number of decimal places.</param>
///// <param name="b">The second value.</param> public static bool AlmostEqualRelative<T>(this Matrix<T> a, Matrix<T> b, int decimalPlaces)
///// <param name="decimalPlaces">The number of decimal places.</param> where T : struct, IEquatable<T>, IFormattable
//public static bool AlmostEqualRelative<T>(this Matrix<T> a, Matrix<T> b, int decimalPlaces) {
// where T : struct, IEquatable<T>, IFormattable return AlmostEqualNormRelative(a.L2Norm(), b.L2Norm(), (a - b).L2Norm(), decimalPlaces);
//{ }
// return AlmostEqualNormRelative(a.L2Norm(), b.L2Norm(), (a - b).L2Norm(), decimalPlaces);
//}
} }
} }

25
src/UnitTests/LinearAlgebraTests/Complex32/MatrixTests.cs

@ -92,13 +92,13 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32
public virtual void CanComputeFrobeniusNorm() public virtual void CanComputeFrobeniusNorm()
{ {
var matrix = TestMatrices["Square3x3"]; var matrix = TestMatrices["Square3x3"];
AssertHelpers.AlmostEqualRelative(11.1427106217473f, matrix.FrobeniusNorm().Real, 6); AssertHelpers.AlmostEqualRelative(11.1427106217473f, matrix.FrobeniusNorm(), 6);
matrix = TestMatrices["Wide2x3"]; matrix = TestMatrices["Wide2x3"];
AssertHelpers.AlmostEqualRelative(5.29055762656452f, matrix.FrobeniusNorm().Real, 6); AssertHelpers.AlmostEqualRelative(5.29055762656452f, matrix.FrobeniusNorm(), 6);
matrix = TestMatrices["Tall3x2"]; matrix = TestMatrices["Tall3x2"];
AssertHelpers.AlmostEqualRelative(7.86574853399217, matrix.FrobeniusNorm().Real, 6); AssertHelpers.AlmostEqualRelative(7.86574853399217, matrix.FrobeniusNorm(), 6);
} }
/// <summary> /// <summary>
@ -108,13 +108,13 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32
public virtual void CanComputeInfinityNorm() public virtual void CanComputeInfinityNorm()
{ {
var matrix = TestMatrices["Square3x3"]; var matrix = TestMatrices["Square3x3"];
AssertHelpers.AlmostEqualRelative(16.7777033f, matrix.InfinityNorm().Real, 6); AssertHelpers.AlmostEqualRelative(16.7777033f, matrix.InfinityNorm(), 6);
matrix = TestMatrices["Wide2x3"]; matrix = TestMatrices["Wide2x3"];
AssertHelpers.AlmostEqualRelative(7.3514039f, matrix.InfinityNorm().Real, 6); AssertHelpers.AlmostEqualRelative(7.3514039f, matrix.InfinityNorm(), 6);
matrix = TestMatrices["Tall3x2"]; matrix = TestMatrices["Tall3x2"];
AssertHelpers.AlmostEqualRelative(10.1023756f, matrix.InfinityNorm().Real, 6); AssertHelpers.AlmostEqualRelative(10.1023756f, matrix.InfinityNorm(), 6);
} }
/// <summary> /// <summary>
@ -124,13 +124,13 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32
public virtual void CanComputeL1Norm() public virtual void CanComputeL1Norm()
{ {
var matrix = TestMatrices["Square3x3"]; var matrix = TestMatrices["Square3x3"];
AssertHelpers.AlmostEqualRelative(12.5401248f, matrix.L1Norm().Real, 7); AssertHelpers.AlmostEqualRelative(12.5401248f, matrix.L1Norm(), 7);
matrix = TestMatrices["Wide2x3"]; matrix = TestMatrices["Wide2x3"];
AssertHelpers.AlmostEqualRelative(5.8647971f, matrix.L1Norm().Real, 7); AssertHelpers.AlmostEqualRelative(5.8647971f, matrix.L1Norm(), 7);
matrix = TestMatrices["Tall3x2"]; matrix = TestMatrices["Tall3x2"];
AssertHelpers.AlmostEqualRelative(9.4933860f, matrix.L1Norm().Real, 7); AssertHelpers.AlmostEqualRelative(9.4933860f, matrix.L1Norm(), 7);
} }
/// <summary> /// <summary>
@ -140,12 +140,13 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32
public virtual void CanComputeL2Norm() public virtual void CanComputeL2Norm()
{ {
var matrix = TestMatrices["Square3x3"]; var matrix = TestMatrices["Square3x3"];
AssertHelpers.AlmostEqualRelative(10.6381752f, matrix.L2Norm().Real, 6); AssertHelpers.AlmostEqualRelative(10.6381752f, matrix.L2Norm(), 6);
matrix = TestMatrices["Wide2x3"]; matrix = TestMatrices["Wide2x3"];
AssertHelpers.AlmostEqualRelative(5.2058554f, matrix.L2Norm().Real, 6); AssertHelpers.AlmostEqualRelative(5.2058554f, matrix.L2Norm(), 6);
matrix = TestMatrices["Tall3x2"]; matrix = TestMatrices["Tall3x2"];
AssertHelpers.AlmostEqualRelative(7.3582664f, matrix.L2Norm().Real, 6); AssertHelpers.AlmostEqualRelative(7.3582664f, matrix.L2Norm(), 6);
} }
/// <summary> /// <summary>

Loading…
Cancel
Save