From 462f0e43a361f96a7e3a6fb54f384ea90fcd4f92 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Thu, 29 May 2014 09:55:20 +0200 Subject: [PATCH] LA: optional argument on EVD decomposition if matrix is known to be symmetric #209 --- .../LinearAlgebra/Complex/DenseMatrix.cs | 4 ++-- .../Complex/Factorization/DenseEvd.cs | 21 ++++++++++------- .../Complex/Factorization/UserEvd.cs | 21 ++++++++++------- src/Numerics/LinearAlgebra/Complex/Matrix.cs | 4 ++-- .../LinearAlgebra/Complex32/DenseMatrix.cs | 4 ++-- .../Complex32/Factorization/DenseEvd.cs | 21 ++++++++++------- .../Complex32/Factorization/UserEvd.cs | 21 ++++++++++------- .../LinearAlgebra/Complex32/Matrix.cs | 4 ++-- .../LinearAlgebra/Double/DenseMatrix.cs | 4 ++-- .../Double/Factorization/DenseEvd.cs | 22 +++++++++++------- .../Double/Factorization/UserEvd.cs | 22 +++++++++++------- src/Numerics/LinearAlgebra/Double/Matrix.cs | 4 ++-- src/Numerics/LinearAlgebra/Matrix.Solve.cs | 2 +- src/Numerics/LinearAlgebra/Options.cs | 23 +++++++++++++++++++ .../LinearAlgebra/Single/DenseMatrix.cs | 4 ++-- .../Single/Factorization/DenseEvd.cs | 22 +++++++++++------- .../Single/Factorization/UserEvd.cs | 22 +++++++++++------- src/Numerics/LinearAlgebra/Single/Matrix.cs | 4 ++-- 18 files changed, 148 insertions(+), 81 deletions(-) diff --git a/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs index 6f86670d..b8e7afc5 100644 --- a/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs @@ -1303,9 +1303,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex return DenseSvd.Create(this, computeVectors); } - public override Evd Evd() + public override Evd Evd(Symmetricity symmetricity = Symmetricity.Unknown) { - return DenseEvd.Create(this); + return DenseEvd.Create(this, symmetricity); } } } diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseEvd.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseEvd.cs index 1f393f50..a8869a5d 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseEvd.cs +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseEvd.cs @@ -62,9 +62,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// the eigenvalue decomposition when the constructor is called and cache it's decomposition. /// /// The matrix to factor. + /// If it is known whether the matrix is symmetric or not the routine can skip checking it itself. /// If is null. /// If EVD algorithm failed to converge with matrix . - public static DenseEvd Create(DenseMatrix matrix) + public static DenseEvd Create(DenseMatrix matrix, Symmetricity symmetricity) { if (matrix.RowCount != matrix.ColumnCount) { @@ -78,14 +79,18 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization var blockDiagonal = new DenseMatrix(order); var eigenValues = new DenseVector(order); - var isSymmetric = true; - - for (var i = 0; isSymmetric && i < order; i++) + bool isSymmetric; + switch (symmetricity) { - for (var j = 0; isSymmetric && j < order; j++) - { - isSymmetric &= matrix.At(i, j) == matrix.At(j, i).Conjugate(); - } + case Symmetricity.ConjugateSymmetric: + isSymmetric = true; + break; + case Symmetricity.Asymmetric: + isSymmetric = false; + break; + default: + isSymmetric = matrix.IsConjugateSymmetric(); + break; } Control.LinearAlgebraProvider.EigenDecomp(isSymmetric, order, matrix.Values, eigenVectors.Values, eigenValues.Values, blockDiagonal.Values); diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/UserEvd.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/UserEvd.cs index 3a145ae3..5b3c2a8b 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/UserEvd.cs +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/UserEvd.cs @@ -62,9 +62,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// the eigenvalue decomposition when the constructor is called and cache it's decomposition. /// /// The matrix to factor. + /// If it is known whether the matrix is symmetric or not the routine can skip checking it itself. /// If is null. /// If EVD algorithm failed to converge with matrix . - public static UserEvd Create(Matrix matrix) + public static UserEvd Create(Matrix matrix, Symmetricity symmetricity) { if (matrix.RowCount != matrix.ColumnCount) { @@ -78,14 +79,18 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization var blockDiagonal = Matrix.Build.SameAs(matrix, order, order); var eigenValues = new DenseVector(order); - var isSymmetric = true; - - for (var i = 0; isSymmetric && i < order; i++) + bool isSymmetric; + switch (symmetricity) { - for (var j = 0; isSymmetric && j < order; j++) - { - isSymmetric &= matrix.At(i, j) == matrix.At(j, i).Conjugate(); - } + case Symmetricity.ConjugateSymmetric: + isSymmetric = true; + break; + case Symmetricity.Asymmetric: + isSymmetric = false; + break; + default: + isSymmetric = matrix.IsConjugateSymmetric(); + break; } if (isSymmetric) diff --git a/src/Numerics/LinearAlgebra/Complex/Matrix.cs b/src/Numerics/LinearAlgebra/Complex/Matrix.cs index e86f28db..83695eb4 100644 --- a/src/Numerics/LinearAlgebra/Complex/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Complex/Matrix.cs @@ -710,9 +710,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex return UserSvd.Create(this, computeVectors); } - public override Evd Evd() + public override Evd Evd(Symmetricity symmetricity = Symmetricity.Unknown) { - return UserEvd.Create(this); + return UserEvd.Create(this, symmetricity); } } } diff --git a/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs index cb261bd4..a756020e 100644 --- a/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs @@ -1300,9 +1300,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 return DenseSvd.Create(this, computeVectors); } - public override Evd Evd() + public override Evd Evd(Symmetricity symmetricity = Symmetricity.Unknown) { - return DenseEvd.Create(this); + return DenseEvd.Create(this, symmetricity); } } } diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseEvd.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseEvd.cs index 1f3fd3a3..b5110bba 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseEvd.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseEvd.cs @@ -63,9 +63,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// the eigenvalue decomposition when the constructor is called and cache it's decomposition. /// /// The matrix to factor. + /// If it is known whether the matrix is symmetric or not the routine can skip checking it itself. /// If is null. /// If EVD algorithm failed to converge with matrix . - public static DenseEvd Create(DenseMatrix matrix) + public static DenseEvd Create(DenseMatrix matrix, Symmetricity symmetricity) { if (matrix.RowCount != matrix.ColumnCount) { @@ -79,14 +80,18 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization var blockDiagonal = new DenseMatrix(order); var eigenValues = new LinearAlgebra.Complex.DenseVector(order); - var isSymmetric = true; - - for (var i = 0; isSymmetric && i < order; i++) + bool isSymmetric; + switch (symmetricity) { - for (var j = 0; isSymmetric && j < order; j++) - { - isSymmetric &= matrix.At(i, j) == matrix.At(j, i).Conjugate(); - } + case Symmetricity.ConjugateSymmetric: + isSymmetric = true; + break; + case Symmetricity.Asymmetric: + isSymmetric = false; + break; + default: + isSymmetric = matrix.IsConjugateSymmetric(); + break; } Control.LinearAlgebraProvider.EigenDecomp(isSymmetric, order, matrix.Values, eigenVectors.Values, eigenValues.Values, blockDiagonal.Values); diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/UserEvd.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/UserEvd.cs index 69078437..88a81ba1 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/UserEvd.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/UserEvd.cs @@ -61,9 +61,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// the eigenvalue decomposition when the constructor is called and cache it's decomposition. /// /// The matrix to factor. + /// If it is known whether the matrix is symmetric or not the routine can skip checking it itself. /// If is null. /// If EVD algorithm failed to converge with matrix . - public static UserEvd Create(Matrix matrix) + public static UserEvd Create(Matrix matrix, Symmetricity symmetricity) { if (matrix.RowCount != matrix.ColumnCount) { @@ -77,14 +78,18 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization var blockDiagonal = Matrix.Build.SameAs(matrix, order, order); var eigenValues = new LinearAlgebra.Complex.DenseVector(order); - var isSymmetric = true; - - for (var i = 0; isSymmetric && i < order; i++) + bool isSymmetric; + switch (symmetricity) { - for (var j = 0; isSymmetric && j < order; j++) - { - isSymmetric &= matrix.At(i, j) == matrix.At(j, i).Conjugate(); - } + case Symmetricity.ConjugateSymmetric: + isSymmetric = true; + break; + case Symmetricity.Asymmetric: + isSymmetric = false; + break; + default: + isSymmetric = matrix.IsConjugateSymmetric(); + break; } if (isSymmetric) diff --git a/src/Numerics/LinearAlgebra/Complex32/Matrix.cs b/src/Numerics/LinearAlgebra/Complex32/Matrix.cs index 4119d872..e8dded99 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Matrix.cs @@ -704,9 +704,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 return UserSvd.Create(this, computeVectors); } - public override Evd Evd() + public override Evd Evd(Symmetricity symmetricity = Symmetricity.Unknown) { - return UserEvd.Create(this); + return UserEvd.Create(this, symmetricity); } } } diff --git a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs index bb9e12d2..bcf9683e 100644 --- a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs @@ -1230,9 +1230,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double return DenseSvd.Create(this, computeVectors); } - public override Evd Evd() + public override Evd Evd(Symmetricity symmetricity = Symmetricity.Unknown) { - return DenseEvd.Create(this); + return DenseEvd.Create(this, symmetricity); } } } diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/DenseEvd.cs b/src/Numerics/LinearAlgebra/Double/Factorization/DenseEvd.cs index c81a4ab9..1f41f894 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/DenseEvd.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/DenseEvd.cs @@ -62,9 +62,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// the eigenvalue decomposition when the constructor is called and cache it's decomposition. /// /// The matrix to factor. + /// If it is known whether the matrix is symmetric or not the routine can skip checking it itself. /// If is null. /// If EVD algorithm failed to converge with matrix . - public static DenseEvd Create(DenseMatrix matrix) + public static DenseEvd Create(DenseMatrix matrix, Symmetricity symmetricity) { if (matrix.RowCount != matrix.ColumnCount) { @@ -78,14 +79,19 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization var blockDiagonal = new DenseMatrix(order); var eigenValues = new LinearAlgebra.Complex.DenseVector(order); - var isSymmetric = true; - - for (var i = 0; isSymmetric && i < order; i++) + bool isSymmetric; + switch (symmetricity) { - for (var j = 0; isSymmetric && j < order; j++) - { - isSymmetric &= matrix.At(i, j) == matrix.At(j, i); - } + case Symmetricity.Symmetric: + case Symmetricity.ConjugateSymmetric: + isSymmetric = true; + break; + case Symmetricity.Asymmetric: + isSymmetric = false; + break; + default: + isSymmetric = matrix.IsSymmetric(); + break; } Control.LinearAlgebraProvider.EigenDecomp(isSymmetric, order, matrix.Values, eigenVectors.Values, eigenValues.Values, blockDiagonal.Values); diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/UserEvd.cs b/src/Numerics/LinearAlgebra/Double/Factorization/UserEvd.cs index 4ee3a7dd..99ac3278 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/UserEvd.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/UserEvd.cs @@ -62,9 +62,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// the eigenvalue decomposition when the constructor is called and cache it's decomposition. /// /// The matrix to factor. + /// If it is known whether the matrix is symmetric or not the routine can skip checking it itself. /// If is null. /// If EVD algorithm failed to converge with matrix . - public static UserEvd Create(Matrix matrix) + public static UserEvd Create(Matrix matrix, Symmetricity symmetricity) { if (matrix.RowCount != matrix.ColumnCount) { @@ -78,14 +79,19 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization var blockDiagonal = Matrix.Build.SameAs(matrix, order, order); var eigenValues = new LinearAlgebra.Complex.DenseVector(order); - var isSymmetric = true; - - for (var i = 0; isSymmetric && i < order; i++) + bool isSymmetric; + switch (symmetricity) { - for (var j = 0; isSymmetric && j < order; j++) - { - isSymmetric &= matrix.At(i, j) == matrix.At(j, i); - } + case Symmetricity.Symmetric: + case Symmetricity.ConjugateSymmetric: + isSymmetric = true; + break; + case Symmetricity.Asymmetric: + isSymmetric = false; + break; + default: + isSymmetric = matrix.IsSymmetric(); + break; } var d = new double[order]; diff --git a/src/Numerics/LinearAlgebra/Double/Matrix.cs b/src/Numerics/LinearAlgebra/Double/Matrix.cs index 7ba8f7e6..086d097a 100644 --- a/src/Numerics/LinearAlgebra/Double/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Double/Matrix.cs @@ -689,9 +689,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double return UserSvd.Create(this, computeVectors); } - public override Evd Evd() + public override Evd Evd(Symmetricity symmetricity = Symmetricity.Unknown) { - return UserEvd.Create(this); + return UserEvd.Create(this, symmetricity); } } } diff --git a/src/Numerics/LinearAlgebra/Matrix.Solve.cs b/src/Numerics/LinearAlgebra/Matrix.Solve.cs index d5298e73..52217fc6 100644 --- a/src/Numerics/LinearAlgebra/Matrix.Solve.cs +++ b/src/Numerics/LinearAlgebra/Matrix.Solve.cs @@ -78,7 +78,7 @@ namespace MathNet.Numerics.LinearAlgebra /// Computes the EVD decomposition for a matrix. /// /// The EVD decomposition object. - public abstract Evd Evd(); + public abstract Evd Evd(Symmetricity symmetricity = Symmetricity.Unknown); diff --git a/src/Numerics/LinearAlgebra/Options.cs b/src/Numerics/LinearAlgebra/Options.cs index d64aac76..a157d41d 100644 --- a/src/Numerics/LinearAlgebra/Options.cs +++ b/src/Numerics/LinearAlgebra/Options.cs @@ -58,4 +58,27 @@ namespace MathNet.Numerics.LinearAlgebra /// Include = 1 } + + public enum Symmetricity + { + /// + /// It is not known yet whether a matrix is symmetric or not. + /// + Unknown = 0, + + /// + /// A matrix is symmetric + /// + Symmetric = 1, + + /// + /// A matrix is complex conjugate symmetric. + /// + ConjugateSymmetric = 2, + + /// + /// A matrix is not symmetric + /// + Asymmetric = 3 + } } diff --git a/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs index 6d61b1e9..5658ad2f 100644 --- a/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs @@ -1230,9 +1230,9 @@ namespace MathNet.Numerics.LinearAlgebra.Single return DenseSvd.Create(this, computeVectors); } - public override Evd Evd() + public override Evd Evd(Symmetricity symmetricity = Symmetricity.Unknown) { - return DenseEvd.Create(this); + return DenseEvd.Create(this, symmetricity); } } } diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/DenseEvd.cs b/src/Numerics/LinearAlgebra/Single/Factorization/DenseEvd.cs index 40b43784..68eec0a4 100644 --- a/src/Numerics/LinearAlgebra/Single/Factorization/DenseEvd.cs +++ b/src/Numerics/LinearAlgebra/Single/Factorization/DenseEvd.cs @@ -62,9 +62,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// the eigenvalue decomposition when the constructor is called and cache it's decomposition. /// /// The matrix to factor. + /// If it is known whether the matrix is symmetric or not the routine can skip checking it itself. /// If is null. /// If EVD algorithm failed to converge with matrix . - public static DenseEvd Create(DenseMatrix matrix) + public static DenseEvd Create(DenseMatrix matrix, Symmetricity symmetricity) { if (matrix.RowCount != matrix.ColumnCount) { @@ -78,14 +79,19 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization var blockDiagonal = new DenseMatrix(order); var eigenValues = new LinearAlgebra.Complex.DenseVector(order); - var isSymmetric = true; - - for (var i = 0; isSymmetric && i < order; i++) + bool isSymmetric; + switch (symmetricity) { - for (var j = 0; isSymmetric && j < order; j++) - { - isSymmetric &= matrix.At(i, j) == matrix.At(j, i); - } + case Symmetricity.Symmetric: + case Symmetricity.ConjugateSymmetric: + isSymmetric = true; + break; + case Symmetricity.Asymmetric: + isSymmetric = false; + break; + default: + isSymmetric = matrix.IsSymmetric(); + break; } Control.LinearAlgebraProvider.EigenDecomp(isSymmetric, order, matrix.Values, eigenVectors.Values, eigenValues.Values, blockDiagonal.Values); diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/UserEvd.cs b/src/Numerics/LinearAlgebra/Single/Factorization/UserEvd.cs index 3d5f10cb..fc06c2e3 100644 --- a/src/Numerics/LinearAlgebra/Single/Factorization/UserEvd.cs +++ b/src/Numerics/LinearAlgebra/Single/Factorization/UserEvd.cs @@ -61,9 +61,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// the eigenvalue decomposition when the constructor is called and cache it's decomposition. /// /// The matrix to factor. + /// If it is known whether the matrix is symmetric or not the routine can skip checking it itself. /// If is null. /// If EVD algorithm failed to converge with matrix . - public static UserEvd Create(Matrix matrix) + public static UserEvd Create(Matrix matrix, Symmetricity symmetricity) { if (matrix.RowCount != matrix.ColumnCount) { @@ -77,14 +78,19 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization var blockDiagonal = Matrix.Build.SameAs(matrix, order, order); var eigenValues = new LinearAlgebra.Complex.DenseVector(order); - var isSymmetric = true; - - for (var i = 0; isSymmetric && i < order; i++) + bool isSymmetric; + switch (symmetricity) { - for (var j = 0; isSymmetric && j < order; j++) - { - isSymmetric &= matrix.At(i, j) == matrix.At(j, i); - } + case Symmetricity.Symmetric: + case Symmetricity.ConjugateSymmetric: + isSymmetric = true; + break; + case Symmetricity.Asymmetric: + isSymmetric = false; + break; + default: + isSymmetric = matrix.IsSymmetric(); + break; } var d = new float[order]; diff --git a/src/Numerics/LinearAlgebra/Single/Matrix.cs b/src/Numerics/LinearAlgebra/Single/Matrix.cs index 349ad3c9..90e8d408 100644 --- a/src/Numerics/LinearAlgebra/Single/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Single/Matrix.cs @@ -689,9 +689,9 @@ namespace MathNet.Numerics.LinearAlgebra.Single return UserSvd.Create(this, computeVectors); } - public override Evd Evd() + public override Evd Evd(Symmetricity symmetricity = Symmetricity.Unknown) { - return UserEvd.Create(this); + return UserEvd.Create(this, symmetricity); } } }