From 9c104be751c260e1f16f89d694a6d038781cbb3c Mon Sep 17 00:00:00 2001
From: Christoph Ruegg <git@cdrnet.ch>
Date: Wed, 25 Sep 2013 22:10:26 +0200
Subject: [PATCH] Basic regression: multiple-response overloads (matrix instead
 of vector)

---
 .../LinearRegression/MultipleRegression.cs    | 36 +++++++++++++++++++
 .../LinearRegression/WeightedRegression.cs    | 24 +++++++++++++
 2 files changed, 60 insertions(+)

diff --git a/src/Numerics/LinearRegression/MultipleRegression.cs b/src/Numerics/LinearRegression/MultipleRegression.cs
index 40655644..e84ed2aa 100644
--- a/src/Numerics/LinearRegression/MultipleRegression.cs
+++ b/src/Numerics/LinearRegression/MultipleRegression.cs
@@ -48,6 +48,18 @@ namespace MathNet.Numerics.LinearRegression
             return x.TransposeThisAndMultiply(x).Cholesky().Solve(x.Transpose()*y);
         }
 
+        /// <summary>
+        /// Find the model parameters β such that X*β with predictor X becomes as close to response Y as possible, with least squares residuals.
+        /// Uses the cholesky decomposition of the normal equations.
+        /// </summary>
+        /// <param name="x">Predictor matrix X</param>
+        /// <param name="y">Response vector Y</param>
+        /// <returns>Best fitting vector for model parameters β</returns>
+        public static Matrix<T> NormalEquations<T>(Matrix<T> x, Matrix<T> y) where T : struct, IEquatable<T>, IFormattable
+        {
+            return x.TransposeThisAndMultiply(x).Cholesky().Solve(x.Transpose() * y);
+        }
+
         /// <summary>
         /// Find the model parameters β such that their linear combination with all predictor-arrays in X become as close to their response in Y as possible, with least squares residuals.
         /// Uses the cholesky decomposition of the normal equations.
@@ -92,6 +104,18 @@ namespace MathNet.Numerics.LinearRegression
             return x.QR().Solve(y);
         }
 
+        /// <summary>
+        /// Find the model parameters β such that X*β with predictor X becomes as close to response Y as possible, with least squares residuals.
+        /// Uses an orthogonal decomposition and is therefore more numerically stable than the normal equations but also slower.
+        /// </summary>
+        /// <param name="x">Predictor matrix X</param>
+        /// <param name="y">Response vector Y</param>
+        /// <returns>Best fitting vector for model parameters β</returns>
+        public static Matrix<T> QR<T>(Matrix<T> x, Matrix<T> y) where T : struct, IEquatable<T>, IFormattable
+        {
+            return x.QR().Solve(y);
+        }
+
         /// <summary>
         /// Find the model parameters β such that their linear combination with all predictor-arrays in X become as close to their response in Y as possible, with least squares residuals.
         /// Uses an orthogonal decomposition and is therefore more numerically stable than the normal equations but also slower.
@@ -135,6 +159,18 @@ namespace MathNet.Numerics.LinearRegression
             return x.Svd().Solve(y);
         }
 
+        /// <summary>
+        /// Find the model parameters β such that X*β with predictor X becomes as close to response Y as possible, with least squares residuals.
+        /// Uses a singular value decomposition and is therefore more numerically stable (especially if ill-conditioned) than the normal equations or QR but also slower.
+        /// </summary>
+        /// <param name="x">Predictor matrix X</param>
+        /// <param name="y">Response vector Y</param>
+        /// <returns>Best fitting vector for model parameters β</returns>
+        public static Matrix<T> Svd<T>(Matrix<T> x, Matrix<T> y) where T : struct, IEquatable<T>, IFormattable
+        {
+            return x.Svd().Solve(y);
+        }
+
         /// <summary>
         /// Find the model parameters β such that their linear combination with all predictor-arrays in X become as close to their response in Y as possible, with least squares residuals.
         /// Uses a singular value decomposition and is therefore more numerically stable (especially if ill-conditioned) than the normal equations or QR but also slower.
diff --git a/src/Numerics/LinearRegression/WeightedRegression.cs b/src/Numerics/LinearRegression/WeightedRegression.cs
index d615ad21..f9aafa1f 100644
--- a/src/Numerics/LinearRegression/WeightedRegression.cs
+++ b/src/Numerics/LinearRegression/WeightedRegression.cs
@@ -45,6 +45,14 @@ namespace MathNet.Numerics.LinearRegression
             return x.TransposeThisAndMultiply(w * x).Cholesky().Solve(x.Transpose() * (w * y));
         }
 
+        /// <summary>
+        /// Weighted Linear Regression using normal equations.
+        /// </summary>
+        public static Matrix<T> Weighted<T>(Matrix<T> x, Matrix<T> y, Matrix<T> w) where T : struct, IEquatable<T>, IFormattable
+        {
+            return x.TransposeThisAndMultiply(w * x).Cholesky().Solve(x.Transpose() * (w * y));
+        }
+
         /// <summary>
         /// Weighted Linear Regression using normal equations.
         /// </summary>
@@ -87,6 +95,22 @@ namespace MathNet.Numerics.LinearRegression
             return Weighted(x, y, w);
         }
 
+        /// <summary>
+        /// Locally-Weighted Linear Regression using normal equations.
+        /// </summary>
+        public static Matrix<T> Local<T>(Matrix<T> x, Matrix<T> y, Vector<T> t, Func<Vector<T>, Vector<T>, T> kernel) where T : struct, IEquatable<T>, IFormattable
+        {
+            // TODO: Kernel definition is a bit weird as it includes computing the difference norm
+            // We can make this more common once we change the norm to always be of type double around LA.
+
+            var w = Matrix<T>.Build.DenseMatrix(x.RowCount, x.RowCount);
+            for (int i = 0; i < x.RowCount; i++)
+            {
+                w.At(i, i, kernel(t, x.Row(i)));
+            }
+            return Weighted(x, y, w);
+        }
+
         public static Func<Vector<double>, Vector<double>, double> GaussianKernel(double radius)
         {
             // TODO: see above...