LA: matrix column/row norms, normalization

src/Numerics/LinearAlgebra/Complex/Matrix.cs

@@ -106,6 +106,104 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             return Math.Sqrt(norm);
         }
 
+        /// <summary>
+        /// Calculates the p-norms of all row vectors.
+        /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm)
+        /// </summary>
+        public override Vector<double> RowNorms(double norm)
+        {
+            if (norm <= 0.0)
+            {
+                throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive);
+            }
+
+            var ret = Vector<double>.Build.Dense(RowCount);
+            if (norm == 2.0)
+            {
+                Storage.FoldRowsUnchecked(ret.Storage, (s, x) => s + x.MagnitudeSquared(), (x, c) => Math.Sqrt(x), ret.Storage, Zeros.AllowSkip);
+            }
+            else if (norm == 1.0)
+            {
+                Storage.FoldRowsUnchecked(ret.Storage, (s, x) => s + x.Magnitude, (x, c) => x, ret.Storage, Zeros.AllowSkip);
+            }
+            else if (double.IsPositiveInfinity(norm))
+            {
+                Storage.FoldRowsUnchecked(ret.Storage, (s, x) => Math.Max(s, x.Magnitude), (x, c) => x, ret.Storage, Zeros.AllowSkip);
+            }
+            else
+            {
+                double invnorm = 1.0/norm;
+                Storage.FoldRowsUnchecked(ret.Storage, (s, x) => s + Math.Pow(x.Magnitude, norm), (x, c) => Math.Pow(x, invnorm), ret.Storage, Zeros.AllowSkip);
+            }
+            return ret;
+        }
+
+        /// <summary>
+        /// Calculates the p-norms of all column vectors.
+        /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm)
+        /// </summary>
+        public override Vector<double> ColumnNorms(double norm)
+        {
+            if (norm <= 0.0)
+            {
+                throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive);
+            }
+
+            var ret = Vector<double>.Build.Dense(ColumnCount);
+            if (norm == 2.0)
+            {
+                Storage.FoldColumnsUnchecked(ret.Storage, (s, x) => s + x.MagnitudeSquared(), (x, c) => Math.Sqrt(x), ret.Storage, Zeros.AllowSkip);
+            }
+            else if (norm == 1.0)
+            {
+                Storage.FoldColumnsUnchecked(ret.Storage, (s, x) => s + x.Magnitude, (x, c) => x, ret.Storage, Zeros.AllowSkip);
+            }
+            else if (double.IsPositiveInfinity(norm))
+            {
+                Storage.FoldColumnsUnchecked(ret.Storage, (s, x) => Math.Max(s, x.Magnitude), (x, c) => x, ret.Storage, Zeros.AllowSkip);
+            }
+            else
+            {
+                double invnorm = 1.0/norm;
+                Storage.FoldColumnsUnchecked(ret.Storage, (s, x) => s + Math.Pow(x.Magnitude, norm), (x, c) => Math.Pow(x, invnorm), ret.Storage, Zeros.AllowSkip);
+            }
+            return ret;
+        }
+
+        /// <summary>
+        /// Normalizes all row vectors to a unit p-norm.
+        /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm)
+        /// </summary>
+        public override sealed Matrix<Complex> NormalizeRows(double norm)
+        {
+            var norminv = ((DenseVectorStorage<double>)RowNorms(norm).Storage).Data;
+            for (int i = 0; i < norminv.Length; i++)
+            {
+                norminv[i] = norminv[i] == 0d ? 1d : 1d/norminv[i];
+            }
+
+            var result = Build.SameAs(this, RowCount, ColumnCount);
+            Storage.MapIndexedTo(result.Storage, (i, j, x) => norminv[i]*x, Zeros.AllowSkip, ExistingData.AssumeZeros);
+            return result;
+        }
+
+        /// <summary>
+        /// Normalizes all column vectors to a unit p-norm.
+        /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm)
+        /// </summary>
+        public override sealed Matrix<Complex> NormalizeColumns(double norm)
+        {
+            var norminv = ((DenseVectorStorage<double>)ColumnNorms(norm).Storage).Data;
+            for (int i = 0; i < norminv.Length; i++)
+            {
+                norminv[i] = norminv[i] == 0d ? 1d : 1d/norminv[i];
+            }
+
+            var result = Build.SameAs(this, RowCount, ColumnCount);
+            Storage.MapIndexedTo(result.Storage, (i, j, x) => norminv[j]*x, Zeros.AllowSkip, ExistingData.AssumeZeros);
+            return result;
+        }
+
         /// <summary>
         /// Returns the conjugate transpose of this matrix.
         /// </summary>

src/Numerics/LinearAlgebra/Complex32/Matrix.cs

@@ -100,6 +100,104 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             return Math.Sqrt(norm);
         }
 
+        /// <summary>
+        /// Calculates the p-norms of all row vectors.
+        /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm)
+        /// </summary>
+        public override Vector<double> RowNorms(double norm)
+        {
+            if (norm <= 0.0)
+            {
+                throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive);
+            }
+
+            var ret = Vector<double>.Build.Dense(RowCount);
+            if (norm == 2.0)
+            {
+                Storage.FoldRowsUnchecked(ret.Storage, (s, x) => s + x.MagnitudeSquared, (x, c) => Math.Sqrt(x), ret.Storage, Zeros.AllowSkip);
+            }
+            else if (norm == 1.0)
+            {
+                Storage.FoldRowsUnchecked(ret.Storage, (s, x) => s + x.Magnitude, (x, c) => x, ret.Storage, Zeros.AllowSkip);
+            }
+            else if (double.IsPositiveInfinity(norm))
+            {
+                Storage.FoldRowsUnchecked(ret.Storage, (s, x) => Math.Max(s, x.Magnitude), (x, c) => x, ret.Storage, Zeros.AllowSkip);
+            }
+            else
+            {
+                double invnorm = 1.0/norm;
+                Storage.FoldRowsUnchecked(ret.Storage, (s, x) => s + Math.Pow(x.Magnitude, norm), (x, c) => Math.Pow(x, invnorm), ret.Storage, Zeros.AllowSkip);
+            }
+            return ret;
+        }
+
+        /// <summary>
+        /// Calculates the p-norms of all column vectors.
+        /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm)
+        /// </summary>
+        public override Vector<double> ColumnNorms(double norm)
+        {
+            if (norm <= 0.0)
+            {
+                throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive);
+            }
+
+            var ret = Vector<double>.Build.Dense(ColumnCount);
+            if (norm == 2.0)
+            {
+                Storage.FoldColumnsUnchecked(ret.Storage, (s, x) => s + x.MagnitudeSquared, (x, c) => Math.Sqrt(x), ret.Storage, Zeros.AllowSkip);
+            }
+            else if (norm == 1.0)
+            {
+                Storage.FoldColumnsUnchecked(ret.Storage, (s, x) => s + x.Magnitude, (x, c) => x, ret.Storage, Zeros.AllowSkip);
+            }
+            else if (double.IsPositiveInfinity(norm))
+            {
+                Storage.FoldColumnsUnchecked(ret.Storage, (s, x) => Math.Max(s, x.Magnitude), (x, c) => x, ret.Storage, Zeros.AllowSkip);
+            }
+            else
+            {
+                double invnorm = 1.0/norm;
+                Storage.FoldColumnsUnchecked(ret.Storage, (s, x) => s + Math.Pow(x.Magnitude, norm), (x, c) => Math.Pow(x, invnorm), ret.Storage, Zeros.AllowSkip);
+            }
+            return ret;
+        }
+
+        /// <summary>
+        /// Normalizes all row vectors to a unit p-norm.
+        /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm)
+        /// </summary>
+        public override sealed Matrix<Complex32> NormalizeRows(double norm)
+        {
+            var norminv = ((DenseVectorStorage<double>)RowNorms(norm).Storage).Data;
+            for (int i = 0; i < norminv.Length; i++)
+            {
+                norminv[i] = norminv[i] == 0d ? 1d : 1d/norminv[i];
+            }
+
+            var result = Build.SameAs(this, RowCount, ColumnCount);
+            Storage.MapIndexedTo(result.Storage, (i, j, x) => ((float)norminv[i])*x, Zeros.AllowSkip, ExistingData.AssumeZeros);
+            return result;
+        }
+
+        /// <summary>
+        /// Normalizes all column vectors to a unit p-norm.
+        /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm)
+        /// </summary>
+        public override sealed Matrix<Complex32> NormalizeColumns(double norm)
+        {
+            var norminv = ((DenseVectorStorage<double>)ColumnNorms(norm).Storage).Data;
+            for (int i = 0; i < norminv.Length; i++)
+            {
+                norminv[i] = norminv[i] == 0d ? 1d : 1d/norminv[i];
+            }
+
+            var result = Build.SameAs(this, RowCount, ColumnCount);
+            Storage.MapIndexedTo(result.Storage, (i, j, x) => ((float)norminv[j])*x, Zeros.AllowSkip, ExistingData.AssumeZeros);
+            return result;
+        }
+
         /// <summary>
         /// Returns the conjugate transpose of this matrix.
         /// </summary>

src/Numerics/LinearAlgebra/Double/Matrix.cs

@@ -98,6 +98,104 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             return Math.Sqrt(norm);
         }
 
+        /// <summary>
+        /// Calculates the p-norms of all row vectors.
+        /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm)
+        /// </summary>
+        public override Vector<double> RowNorms(double norm)
+        {
+            if (norm <= 0.0)
+            {
+                throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive);
+            }
+
+            var ret = Vector<double>.Build.Dense(RowCount);
+            if (norm == 2.0)
+            {
+                Storage.FoldRowsUnchecked(ret.Storage, (s, x) => s + x*x, (x, c) => Math.Sqrt(x), ret.Storage, Zeros.AllowSkip);
+            }
+            else if (norm == 1.0)
+            {
+                Storage.FoldRowsUnchecked(ret.Storage, (s, x) => s + Math.Abs(x), (x, c) => x, ret.Storage, Zeros.AllowSkip);
+            }
+            else if (double.IsPositiveInfinity(norm))
+            {
+                Storage.FoldRowsUnchecked(ret.Storage, (s, x) => Math.Max(s, Math.Abs(x)), (x, c) => x, ret.Storage, Zeros.AllowSkip);
+            }
+            else
+            {
+                double invnorm = 1.0/norm;
+                Storage.FoldRowsUnchecked(ret.Storage, (s, x) => s + Math.Pow(Math.Abs(x), norm), (x, c) => Math.Pow(x, invnorm), ret.Storage, Zeros.AllowSkip);
+            }
+            return ret;
+        }
+
+        /// <summary>
+        /// Calculates the p-norms of all column vectors.
+        /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm)
+        /// </summary>
+        public override Vector<double> ColumnNorms(double norm)
+        {
+            if (norm <= 0.0)
+            {
+                throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive);
+            }
+
+            var ret = Vector<double>.Build.Dense(ColumnCount);
+            if (norm == 2.0)
+            {
+                Storage.FoldColumnsUnchecked(ret.Storage, (s, x) => s + x*x, (x, c) => Math.Sqrt(x), ret.Storage, Zeros.AllowSkip);
+            }
+            else if (norm == 1.0)
+            {
+                Storage.FoldColumnsUnchecked(ret.Storage, (s, x) => s + Math.Abs(x), (x, c) => x, ret.Storage, Zeros.AllowSkip);
+            }
+            else if (double.IsPositiveInfinity(norm))
+            {
+                Storage.FoldColumnsUnchecked(ret.Storage, (s, x) => Math.Max(s, Math.Abs(x)), (x, c) => x, ret.Storage, Zeros.AllowSkip);
+            }
+            else
+            {
+                double invnorm = 1.0/norm;
+                Storage.FoldColumnsUnchecked(ret.Storage, (s, x) => s + Math.Pow(Math.Abs(x), norm), (x, c) => Math.Pow(x, invnorm), ret.Storage, Zeros.AllowSkip);
+            }
+            return ret;
+        }
+
+        /// <summary>
+        /// Normalizes all row vectors to a unit p-norm.
+        /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm)
+        /// </summary>
+        public override sealed Matrix<double> NormalizeRows(double norm)
+        {
+            var norminv = ((DenseVectorStorage<double>)RowNorms(norm).Storage).Data;
+            for (int i = 0; i < norminv.Length; i++)
+            {
+                norminv[i] = norminv[i] == 0d ? 1d : 1d/norminv[i];
+            }
+
+            var result = Build.SameAs(this, RowCount, ColumnCount);
+            Storage.MapIndexedTo(result.Storage, (i, j, x) => norminv[i]*x, Zeros.AllowSkip, ExistingData.AssumeZeros);
+            return result;
+        }
+
+        /// <summary>
+        /// Normalizes all column vectors to a unit p-norm.
+        /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm)
+        /// </summary>
+        public override sealed Matrix<double> NormalizeColumns(double norm)
+        {
+            var norminv = ((DenseVectorStorage<double>)ColumnNorms(norm).Storage).Data;
+            for (int i = 0; i < norminv.Length; i++)
+            {
+                norminv[i] = norminv[i] == 0d ? 1d : 1d/norminv[i];
+            }
+
+            var result = Build.SameAs(this, RowCount, ColumnCount);
+            Storage.MapIndexedTo(result.Storage, (i, j, x) => norminv[j]*x, Zeros.AllowSkip, ExistingData.AssumeZeros);
+            return result;
+        }
+
         /// <summary>
         /// Returns the conjugate transpose of this matrix.
         /// </summary>

src/Numerics/LinearAlgebra/Matrix.Arithmetic.cs

@@ -1627,50 +1627,6 @@ namespace MathNet.Numerics.LinearAlgebra
             }
         }
 
-        /// <summary>
-        /// Normalizes the columns of a matrix.
-        /// </summary>
-        /// <param name="p">The norm under which to normalize the columns under.</param>
-        /// <returns>A normalized version of the matrix.</returns>
-        /// <exception cref="ArgumentOutOfRangeException">If the parameter p is not positive.</exception>
-        public Matrix<T> NormalizeColumns(int p)
-        {
-            if (p < 1)
-            {
-                throw new ArgumentOutOfRangeException("p", Resources.ArgumentMustBePositive);
-            }
-
-            var result = Build.SameAs(this);
-            for (var index = 0; index < ColumnCount; index++)
-            {
-                result.SetColumn(index, Column(index).Normalize(p));
-            }
-
-            return result;
-        }
-
-        /// <summary>
-        /// Normalizes the rows of a matrix.
-        /// </summary>
-        /// <param name="p">The norm under which to normalize the rows under.</param>
-        /// <returns>A normalized version of the matrix.</returns>
-        /// <exception cref="ArgumentOutOfRangeException">If the parameter p is not positive.</exception>
-        public Matrix<T> NormalizeRows(int p)
-        {
-            if (p < 1)
-            {
-                throw new ArgumentOutOfRangeException("p", Resources.ArgumentMustBePositive);
-            }
-
-            var ret = Build.SameAs(this);
-            for (var index = 0; index < RowCount; index++)
-            {
-                ret.SetRow(index, Row(index).Normalize(p));
-            }
-
-            return ret;
-        }
-
         /// <summary>Calculates the induced L1 norm of this matrix.</summary>
         /// <returns>The maximum absolute column sum of the matrix.</returns>
         public abstract double L1Norm();
@@ -1694,7 +1650,29 @@ namespace MathNet.Numerics.LinearAlgebra
         /// <returns>The square root of the sum of the squared values.</returns>
         public abstract double FrobeniusNorm();
 
+        /// <summary>
+        /// Calculates the p-norms of all row vectors.
+        /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm)
+        /// </summary>
+        public abstract Vector<double> RowNorms(double norm);
 
+        /// <summary>
+        /// Calculates the p-norms of all column vectors.
+        /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm)
+        /// </summary>
+        public abstract Vector<double> ColumnNorms(double norm);
+
+        /// <summary>
+        /// Normalizes all row vectors to a unit p-norm.
+        /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm)
+        /// </summary>
+        public abstract Matrix<T> NormalizeRows(double norm);
+
+        /// <summary>
+        /// Normalizes all column vectors to a unit p-norm.
+        /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm)
+        /// </summary>
+        public abstract Matrix<T> NormalizeColumns(double norm);
 
         #region Exceptions - possibly move elsewhere?
 

src/Numerics/LinearAlgebra/Single/Matrix.cs

@@ -98,6 +98,104 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             return Math.Sqrt(norm);
         }
 
+        /// <summary>
+        /// Calculates the p-norms of all row vectors.
+        /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm)
+        /// </summary>
+        public override Vector<double> RowNorms(double norm)
+        {
+            if (norm <= 0.0)
+            {
+                throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive);
+            }
+
+            var ret = Vector<double>.Build.Dense(RowCount);
+            if (norm == 2.0)
+            {
+                Storage.FoldRowsUnchecked(ret.Storage, (s, x) => s + x*x, (x, c) => Math.Sqrt(x), ret.Storage, Zeros.AllowSkip);
+            }
+            else if (norm == 1.0)
+            {
+                Storage.FoldRowsUnchecked(ret.Storage, (s, x) => s + Math.Abs(x), (x, c) => x, ret.Storage, Zeros.AllowSkip);
+            }
+            else if (double.IsPositiveInfinity(norm))
+            {
+                Storage.FoldRowsUnchecked(ret.Storage, (s, x) => Math.Max(s, Math.Abs(x)), (x, c) => x, ret.Storage, Zeros.AllowSkip);
+            }
+            else
+            {
+                double invnorm = 1.0/norm;
+                Storage.FoldRowsUnchecked(ret.Storage, (s, x) => s + Math.Pow(Math.Abs(x), norm), (x, c) => Math.Pow(x, invnorm), ret.Storage, Zeros.AllowSkip);
+            }
+            return ret;
+        }
+
+        /// <summary>
+        /// Calculates the p-norms of all column vectors.
+        /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm)
+        /// </summary>
+        public override Vector<double> ColumnNorms(double norm)
+        {
+            if (norm <= 0.0)
+            {
+                throw new ArgumentOutOfRangeException("norm", Resources.ArgumentMustBePositive);
+            }
+
+            var ret = Vector<double>.Build.Dense(ColumnCount);
+            if (norm == 2.0)
+            {
+                Storage.FoldColumnsUnchecked(ret.Storage, (s, x) => s + x*x, (x, c) => Math.Sqrt(x), ret.Storage, Zeros.AllowSkip);
+            }
+            else if (norm == 1.0)
+            {
+                Storage.FoldColumnsUnchecked(ret.Storage, (s, x) => s + Math.Abs(x), (x, c) => x, ret.Storage, Zeros.AllowSkip);
+            }
+            else if (double.IsPositiveInfinity(norm))
+            {
+                Storage.FoldColumnsUnchecked(ret.Storage, (s, x) => Math.Max(s, Math.Abs(x)), (x, c) => x, ret.Storage, Zeros.AllowSkip);
+            }
+            else
+            {
+                double invnorm = 1.0/norm;
+                Storage.FoldColumnsUnchecked(ret.Storage, (s, x) => s + Math.Pow(Math.Abs(x), norm), (x, c) => Math.Pow(x, invnorm), ret.Storage, Zeros.AllowSkip);
+            }
+            return ret;
+        }
+
+        /// <summary>
+        /// Normalizes all row vectors to a unit p-norm.
+        /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm)
+        /// </summary>
+        public override sealed Matrix<float> NormalizeRows(double norm)
+        {
+            var norminv = ((DenseVectorStorage<double>)RowNorms(norm).Storage).Data;
+            for (int i = 0; i < norminv.Length; i++)
+            {
+                norminv[i] = norminv[i] == 0d ? 1d : 1d/norminv[i];
+            }
+
+            var result = Build.SameAs(this, RowCount, ColumnCount);
+            Storage.MapIndexedTo(result.Storage, (i, j, x) => (float)norminv[i]*x, Zeros.AllowSkip, ExistingData.AssumeZeros);
+            return result;
+        }
+
+        /// <summary>
+        /// Normalizes all column vectors to a unit p-norm.
+        /// Typical values for p are 1.0 (L1, Manhattan norm), 2.0 (L2, Euclidean norm) and positive infinity (infinity norm)
+        /// </summary>
+        public override sealed Matrix<float> NormalizeColumns(double norm)
+        {
+            var norminv = ((DenseVectorStorage<double>)ColumnNorms(norm).Storage).Data;
+            for (int i = 0; i < norminv.Length; i++)
+            {
+                norminv[i] = norminv[i] == 0d ? 1d : 1d/norminv[i];
+            }
+
+            var result = Build.SameAs(this, RowCount, ColumnCount);
+            Storage.MapIndexedTo(result.Storage, (i, j, x) => (float)norminv[j]*x, Zeros.AllowSkip, ExistingData.AssumeZeros);
+            return result;
+        }
+
         /// <summary>
         /// Returns the conjugate transpose of this matrix.
         /// </summary>