[8323] | 1 | #region License Information
|
---|
| 2 | /* HeuristicLab
|
---|
[11594] | 3 | * Copyright (C) 2002-2014 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
|
---|
[8323] | 4 | *
|
---|
| 5 | * This file is part of HeuristicLab.
|
---|
| 6 | *
|
---|
| 7 | * HeuristicLab is free software: you can redistribute it and/or modify
|
---|
| 8 | * it under the terms of the GNU General Public License as published by
|
---|
| 9 | * the Free Software Foundation, either version 3 of the License, or
|
---|
| 10 | * (at your option) any later version.
|
---|
| 11 | *
|
---|
| 12 | * HeuristicLab is distributed in the hope that it will be useful,
|
---|
| 13 | * but WITHOUT ANY WARRANTY; without even the implied warranty of
|
---|
| 14 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
---|
| 15 | * GNU General Public License for more details.
|
---|
| 16 | *
|
---|
| 17 | * You should have received a copy of the GNU General Public License
|
---|
| 18 | * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
|
---|
| 19 | */
|
---|
| 20 | #endregion
|
---|
| 21 |
|
---|
| 22 | using System;
|
---|
| 23 | using System.Collections.Generic;
|
---|
| 24 | using System.Linq;
|
---|
| 25 | using HeuristicLab.Common;
|
---|
| 26 | using HeuristicLab.Core;
|
---|
| 27 | using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
|
---|
| 28 | using HeuristicLab.Problems.DataAnalysis;
|
---|
| 29 |
|
---|
[8371] | 30 | namespace HeuristicLab.Algorithms.DataAnalysis {
|
---|
[8323] | 31 | /// <summary>
|
---|
| 32 | /// Represents a Gaussian process model.
|
---|
| 33 | /// </summary>
|
---|
| 34 | [StorableClass]
|
---|
| 35 | [Item("GaussianProcessModel", "Represents a Gaussian process posterior.")]
|
---|
| 36 | public sealed class GaussianProcessModel : NamedItem, IGaussianProcessModel {
|
---|
| 37 | [Storable]
|
---|
| 38 | private double negativeLogLikelihood;
|
---|
| 39 | public double NegativeLogLikelihood {
|
---|
| 40 | get { return negativeLogLikelihood; }
|
---|
| 41 | }
|
---|
| 42 |
|
---|
| 43 | [Storable]
|
---|
[8484] | 44 | private double[] hyperparameterGradients;
|
---|
| 45 | public double[] HyperparameterGradients {
|
---|
| 46 | get {
|
---|
| 47 | var copy = new double[hyperparameterGradients.Length];
|
---|
| 48 | Array.Copy(hyperparameterGradients, copy, copy.Length);
|
---|
| 49 | return copy;
|
---|
| 50 | }
|
---|
| 51 | }
|
---|
| 52 |
|
---|
| 53 | [Storable]
|
---|
[8323] | 54 | private ICovarianceFunction covarianceFunction;
|
---|
| 55 | public ICovarianceFunction CovarianceFunction {
|
---|
| 56 | get { return covarianceFunction; }
|
---|
| 57 | }
|
---|
| 58 | [Storable]
|
---|
| 59 | private IMeanFunction meanFunction;
|
---|
| 60 | public IMeanFunction MeanFunction {
|
---|
| 61 | get { return meanFunction; }
|
---|
| 62 | }
|
---|
| 63 | [Storable]
|
---|
| 64 | private string targetVariable;
|
---|
| 65 | public string TargetVariable {
|
---|
| 66 | get { return targetVariable; }
|
---|
| 67 | }
|
---|
| 68 | [Storable]
|
---|
| 69 | private string[] allowedInputVariables;
|
---|
| 70 | public string[] AllowedInputVariables {
|
---|
| 71 | get { return allowedInputVariables; }
|
---|
| 72 | }
|
---|
| 73 |
|
---|
| 74 | [Storable]
|
---|
| 75 | private double[] alpha;
|
---|
| 76 | [Storable]
|
---|
| 77 | private double sqrSigmaNoise;
|
---|
[8582] | 78 | public double SigmaNoise {
|
---|
| 79 | get { return Math.Sqrt(sqrSigmaNoise); }
|
---|
| 80 | }
|
---|
[8323] | 81 |
|
---|
| 82 | [Storable]
|
---|
[8982] | 83 | private double[] meanParameter;
|
---|
| 84 | [Storable]
|
---|
| 85 | private double[] covarianceParameter;
|
---|
| 86 |
|
---|
| 87 | [Storable]
|
---|
[8323] | 88 | private double[,] l;
|
---|
| 89 |
|
---|
| 90 | [Storable]
|
---|
| 91 | private double[,] x;
|
---|
| 92 | [Storable]
|
---|
[8463] | 93 | private Scaling inputScaling;
|
---|
[8323] | 94 |
|
---|
| 95 |
|
---|
| 96 | [StorableConstructor]
|
---|
| 97 | private GaussianProcessModel(bool deserializing) : base(deserializing) { }
|
---|
| 98 | private GaussianProcessModel(GaussianProcessModel original, Cloner cloner)
|
---|
| 99 | : base(original, cloner) {
|
---|
| 100 | this.meanFunction = cloner.Clone(original.meanFunction);
|
---|
| 101 | this.covarianceFunction = cloner.Clone(original.covarianceFunction);
|
---|
[8463] | 102 | this.inputScaling = cloner.Clone(original.inputScaling);
|
---|
[8323] | 103 | this.negativeLogLikelihood = original.negativeLogLikelihood;
|
---|
| 104 | this.targetVariable = original.targetVariable;
|
---|
[8416] | 105 | this.sqrSigmaNoise = original.sqrSigmaNoise;
|
---|
[8982] | 106 | if (original.meanParameter != null) {
|
---|
| 107 | this.meanParameter = (double[])original.meanParameter.Clone();
|
---|
| 108 | }
|
---|
| 109 | if (original.covarianceParameter != null) {
|
---|
| 110 | this.covarianceParameter = (double[])original.covarianceParameter.Clone();
|
---|
| 111 | }
|
---|
[8416] | 112 |
|
---|
| 113 | // shallow copies of arrays because they cannot be modified
|
---|
[8323] | 114 | this.allowedInputVariables = original.allowedInputVariables;
|
---|
| 115 | this.alpha = original.alpha;
|
---|
| 116 | this.l = original.l;
|
---|
| 117 | this.x = original.x;
|
---|
| 118 | }
|
---|
| 119 | public GaussianProcessModel(Dataset ds, string targetVariable, IEnumerable<string> allowedInputVariables, IEnumerable<int> rows,
|
---|
| 120 | IEnumerable<double> hyp, IMeanFunction meanFunction, ICovarianceFunction covarianceFunction)
|
---|
| 121 | : base() {
|
---|
| 122 | this.name = ItemName;
|
---|
| 123 | this.description = ItemDescription;
|
---|
[8416] | 124 | this.meanFunction = (IMeanFunction)meanFunction.Clone();
|
---|
| 125 | this.covarianceFunction = (ICovarianceFunction)covarianceFunction.Clone();
|
---|
[8323] | 126 | this.targetVariable = targetVariable;
|
---|
| 127 | this.allowedInputVariables = allowedInputVariables.ToArray();
|
---|
| 128 |
|
---|
| 129 |
|
---|
[8416] | 130 | int nVariables = this.allowedInputVariables.Length;
|
---|
[8982] | 131 | meanParameter = hyp
|
---|
[8416] | 132 | .Take(this.meanFunction.GetNumberOfParameters(nVariables))
|
---|
[8982] | 133 | .ToArray();
|
---|
| 134 |
|
---|
| 135 | covarianceParameter = hyp.Skip(this.meanFunction.GetNumberOfParameters(nVariables))
|
---|
| 136 | .Take(this.covarianceFunction.GetNumberOfParameters(nVariables))
|
---|
| 137 | .ToArray();
|
---|
[8473] | 138 | sqrSigmaNoise = Math.Exp(2.0 * hyp.Last());
|
---|
[8416] | 139 |
|
---|
| 140 | CalculateModel(ds, rows);
|
---|
[8323] | 141 | }
|
---|
| 142 |
|
---|
[8416] | 143 | private void CalculateModel(Dataset ds, IEnumerable<int> rows) {
|
---|
[8463] | 144 | inputScaling = new Scaling(ds, allowedInputVariables, rows);
|
---|
| 145 | x = AlglibUtil.PrepareAndScaleInputMatrix(ds, allowedInputVariables, rows, inputScaling);
|
---|
[8473] | 146 | var y = ds.GetDoubleValues(targetVariable, rows);
|
---|
[8323] | 147 |
|
---|
| 148 | int n = x.GetLength(0);
|
---|
| 149 | l = new double[n, n];
|
---|
| 150 |
|
---|
| 151 | // calculate means and covariances
|
---|
[8982] | 152 | var mean = meanFunction.GetParameterizedMeanFunction(meanParameter, Enumerable.Range(0, x.GetLength(1)));
|
---|
| 153 | double[] m = Enumerable.Range(0, x.GetLength(0))
|
---|
| 154 | .Select(r => mean.Mean(x, r))
|
---|
| 155 | .ToArray();
|
---|
| 156 |
|
---|
[9357] | 157 | var cov = covarianceFunction.GetParameterizedCovarianceFunction(covarianceParameter, Enumerable.Range(0, x.GetLength(1)));
|
---|
[8323] | 158 | for (int i = 0; i < n; i++) {
|
---|
| 159 | for (int j = i; j < n; j++) {
|
---|
[8982] | 160 | l[j, i] = cov.Covariance(x, i, j) / sqrSigmaNoise;
|
---|
[8323] | 161 | if (j == i) l[j, i] += 1.0;
|
---|
| 162 | }
|
---|
| 163 | }
|
---|
| 164 |
|
---|
[8982] | 165 |
|
---|
[8323] | 166 | // cholesky decomposition
|
---|
| 167 | int info;
|
---|
| 168 | alglib.densesolverreport denseSolveRep;
|
---|
| 169 |
|
---|
| 170 | var res = alglib.trfac.spdmatrixcholesky(ref l, n, false);
|
---|
[8473] | 171 | if (!res) throw new ArgumentException("Matrix is not positive semidefinite");
|
---|
[8323] | 172 |
|
---|
| 173 | // calculate sum of diagonal elements for likelihood
|
---|
| 174 | double diagSum = Enumerable.Range(0, n).Select(i => Math.Log(l[i, i])).Sum();
|
---|
| 175 |
|
---|
| 176 | // solve for alpha
|
---|
| 177 | double[] ym = y.Zip(m, (a, b) => a - b).ToArray();
|
---|
| 178 |
|
---|
| 179 | alglib.spdmatrixcholeskysolve(l, n, false, ym, out info, out denseSolveRep, out alpha);
|
---|
| 180 | for (int i = 0; i < alpha.Length; i++)
|
---|
| 181 | alpha[i] = alpha[i] / sqrSigmaNoise;
|
---|
| 182 | negativeLogLikelihood = 0.5 * Util.ScalarProd(ym, alpha) + diagSum + (n / 2.0) * Math.Log(2.0 * Math.PI * sqrSigmaNoise);
|
---|
| 183 |
|
---|
| 184 | // derivatives
|
---|
| 185 | int nAllowedVariables = x.GetLength(1);
|
---|
| 186 |
|
---|
[8463] | 187 | alglib.matinvreport matInvRep;
|
---|
[8475] | 188 | double[,] lCopy = new double[l.GetLength(0), l.GetLength(1)];
|
---|
| 189 | Array.Copy(l, lCopy, lCopy.Length);
|
---|
[8323] | 190 |
|
---|
[8475] | 191 | alglib.spdmatrixcholeskyinverse(ref lCopy, n, false, out info, out matInvRep);
|
---|
[8463] | 192 | if (info != 1) throw new ArgumentException("Can't invert matrix to calculate gradients.");
|
---|
[8323] | 193 | for (int i = 0; i < n; i++) {
|
---|
[8463] | 194 | for (int j = 0; j <= i; j++)
|
---|
[8475] | 195 | lCopy[i, j] = lCopy[i, j] / sqrSigmaNoise - alpha[i] * alpha[j];
|
---|
[8323] | 196 | }
|
---|
| 197 |
|
---|
[8475] | 198 | double noiseGradient = sqrSigmaNoise * Enumerable.Range(0, n).Select(i => lCopy[i, i]).Sum();
|
---|
[8323] | 199 |
|
---|
| 200 | double[] meanGradients = new double[meanFunction.GetNumberOfParameters(nAllowedVariables)];
|
---|
[8982] | 201 | for (int k = 0; k < meanGradients.Length; k++) {
|
---|
| 202 | var meanGrad = Enumerable.Range(0, alpha.Length)
|
---|
| 203 | .Select(r => mean.Gradient(x, r, k));
|
---|
| 204 | meanGradients[k] = -Util.ScalarProd(meanGrad, alpha);
|
---|
[8323] | 205 | }
|
---|
| 206 |
|
---|
| 207 | double[] covGradients = new double[covarianceFunction.GetNumberOfParameters(nAllowedVariables)];
|
---|
[8366] | 208 | if (covGradients.Length > 0) {
|
---|
| 209 | for (int i = 0; i < n; i++) {
|
---|
[8484] | 210 | for (int j = 0; j < i; j++) {
|
---|
[8982] | 211 | var g = cov.CovarianceGradient(x, i, j).ToArray();
|
---|
[8484] | 212 | for (int k = 0; k < covGradients.Length; k++) {
|
---|
| 213 | covGradients[k] += lCopy[i, j] * g[k];
|
---|
[8366] | 214 | }
|
---|
[8323] | 215 | }
|
---|
[8484] | 216 |
|
---|
[8982] | 217 | var gDiag = cov.CovarianceGradient(x, i, i).ToArray();
|
---|
[8484] | 218 | for (int k = 0; k < covGradients.Length; k++) {
|
---|
| 219 | // diag
|
---|
| 220 | covGradients[k] += 0.5 * lCopy[i, i] * gDiag[k];
|
---|
| 221 | }
|
---|
[8323] | 222 | }
|
---|
| 223 | }
|
---|
| 224 |
|
---|
[8484] | 225 | hyperparameterGradients =
|
---|
[8473] | 226 | meanGradients
|
---|
| 227 | .Concat(covGradients)
|
---|
| 228 | .Concat(new double[] { noiseGradient }).ToArray();
|
---|
[8484] | 229 |
|
---|
[8323] | 230 | }
|
---|
| 231 |
|
---|
| 232 |
|
---|
| 233 | public override IDeepCloneable Clone(Cloner cloner) {
|
---|
| 234 | return new GaussianProcessModel(this, cloner);
|
---|
| 235 | }
|
---|
| 236 |
|
---|
[8982] | 237 | // is called by the solution creator to set all parameter values of the covariance and mean function
|
---|
| 238 | // to the optimized values (necessary to make the values visible in the GUI)
|
---|
| 239 | public void FixParameters() {
|
---|
| 240 | covarianceFunction.SetParameter(covarianceParameter);
|
---|
| 241 | meanFunction.SetParameter(meanParameter);
|
---|
| 242 | covarianceParameter = new double[0];
|
---|
| 243 | meanParameter = new double[0];
|
---|
| 244 | }
|
---|
| 245 |
|
---|
[8323] | 246 | #region IRegressionModel Members
|
---|
| 247 | public IEnumerable<double> GetEstimatedValues(Dataset dataset, IEnumerable<int> rows) {
|
---|
| 248 | return GetEstimatedValuesHelper(dataset, rows);
|
---|
| 249 | }
|
---|
| 250 | public GaussianProcessRegressionSolution CreateRegressionSolution(IRegressionProblemData problemData) {
|
---|
[8528] | 251 | return new GaussianProcessRegressionSolution(this, new RegressionProblemData(problemData));
|
---|
[8323] | 252 | }
|
---|
| 253 | IRegressionSolution IRegressionModel.CreateRegressionSolution(IRegressionProblemData problemData) {
|
---|
| 254 | return CreateRegressionSolution(problemData);
|
---|
| 255 | }
|
---|
| 256 | #endregion
|
---|
| 257 |
|
---|
[8623] | 258 |
|
---|
[8323] | 259 | private IEnumerable<double> GetEstimatedValuesHelper(Dataset dataset, IEnumerable<int> rows) {
|
---|
[8463] | 260 | var newX = AlglibUtil.PrepareAndScaleInputMatrix(dataset, allowedInputVariables, rows, inputScaling);
|
---|
[8323] | 261 | int newN = newX.GetLength(0);
|
---|
| 262 | int n = x.GetLength(0);
|
---|
| 263 | var Ks = new double[newN, n];
|
---|
[8982] | 264 | var mean = meanFunction.GetParameterizedMeanFunction(meanParameter, Enumerable.Range(0, newX.GetLength(1)));
|
---|
| 265 | var ms = Enumerable.Range(0, newX.GetLength(0))
|
---|
| 266 | .Select(r => mean.Mean(newX, r))
|
---|
| 267 | .ToArray();
|
---|
[9358] | 268 | var cov = covarianceFunction.GetParameterizedCovarianceFunction(covarianceParameter, Enumerable.Range(0, newX.GetLength(1)));
|
---|
[8323] | 269 | for (int i = 0; i < newN; i++) {
|
---|
| 270 | for (int j = 0; j < n; j++) {
|
---|
[8982] | 271 | Ks[i, j] = cov.CrossCovariance(x, newX, j, i);
|
---|
[8323] | 272 | }
|
---|
| 273 | }
|
---|
| 274 |
|
---|
[8463] | 275 | return Enumerable.Range(0, newN)
|
---|
[8473] | 276 | .Select(i => ms[i] + Util.ScalarProd(Util.GetRow(Ks, i), alpha));
|
---|
[8323] | 277 | }
|
---|
[8473] | 278 |
|
---|
| 279 | public IEnumerable<double> GetEstimatedVariance(Dataset dataset, IEnumerable<int> rows) {
|
---|
| 280 | var newX = AlglibUtil.PrepareAndScaleInputMatrix(dataset, allowedInputVariables, rows, inputScaling);
|
---|
| 281 | int newN = newX.GetLength(0);
|
---|
| 282 | int n = x.GetLength(0);
|
---|
| 283 |
|
---|
| 284 | var kss = new double[newN];
|
---|
| 285 | double[,] sWKs = new double[n, newN];
|
---|
[9357] | 286 | var cov = covarianceFunction.GetParameterizedCovarianceFunction(covarianceParameter, Enumerable.Range(0, x.GetLength(1)));
|
---|
[8473] | 287 |
|
---|
| 288 | // for stddev
|
---|
| 289 | for (int i = 0; i < newN; i++)
|
---|
[8982] | 290 | kss[i] = cov.Covariance(newX, i, i);
|
---|
[8473] | 291 |
|
---|
[8475] | 292 | for (int i = 0; i < newN; i++) {
|
---|
| 293 | for (int j = 0; j < n; j++) {
|
---|
[8982] | 294 | sWKs[j, i] = cov.CrossCovariance(x, newX, j, i) / Math.Sqrt(sqrSigmaNoise);
|
---|
[8473] | 295 | }
|
---|
| 296 | }
|
---|
| 297 |
|
---|
| 298 | // for stddev
|
---|
[8484] | 299 | alglib.ablas.rmatrixlefttrsm(n, newN, l, 0, 0, false, false, 0, ref sWKs, 0, 0);
|
---|
[8473] | 300 |
|
---|
| 301 | for (int i = 0; i < newN; i++) {
|
---|
[8484] | 302 | var sumV = Util.ScalarProd(Util.GetCol(sWKs, i), Util.GetCol(sWKs, i));
|
---|
[8475] | 303 | kss[i] -= sumV;
|
---|
| 304 | if (kss[i] < 0) kss[i] = 0;
|
---|
[8473] | 305 | }
|
---|
[8475] | 306 | return kss;
|
---|
[8473] | 307 | }
|
---|
[8323] | 308 | }
|
---|
| 309 | }
|
---|