Free cookie consent management tool by TermsFeed Policy Generator

source: branches/HeuristicLab.Problems.GaussianProcessTuning/ILNumerics.2.14.4735.573/Algorithms/MachineLearning/ILKRR.cs @ 10617

Last change on this file since 10617 was 9102, checked in by gkronber, 12 years ago

#1967: ILNumerics source for experimentation
File size: 12.7 KB
Line 
1///
2///    This file is part of ILNumerics Community Edition.
3///
4///    ILNumerics Community Edition - high performance computing for applications.
5///    Copyright (C) 2006 - 2012 Haymo Kutschbach, http://ilnumerics.net
6///
7///    ILNumerics Community Edition is free software: you can redistribute it and/or modify
8///    it under the terms of the GNU General Public License version 3 as published by
9///    the Free Software Foundation.
10///
11///    ILNumerics Community Edition is distributed in the hope that it will be useful,
12///    but WITHOUT ANY WARRANTY; without even the implied warranty of
13///    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14///    GNU General Public License for more details.
15///
16///    You should have received a copy of the GNU General Public License
17///    along with ILNumerics Community Edition. See the file License.txt in the root
18///    of your distribution package. If not, see <http://www.gnu.org/licenses/>.
19///
20///    In addition this software uses the following components and/or licenses:
21///
22///    =================================================================================
23///    The Open Toolkit Library License
24///   
25///    Copyright (c) 2006 - 2009 the Open Toolkit library.
26///   
27///    Permission is hereby granted, free of charge, to any person obtaining a copy
28///    of this software and associated documentation files (the "Software"), to deal
29///    in the Software without restriction, including without limitation the rights to
30///    use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of
31///    the Software, and to permit persons to whom the Software is furnished to do
32///    so, subject to the following conditions:
33///
34///    The above copyright notice and this permission notice shall be included in all
35///    copies or substantial portions of the Software.
36///
37///    =================================================================================
38///   
39
40using System;
41using System.Collections.Generic;
42using System.Linq;
43using System.Text;
44
45namespace ILNumerics {
46    public partial class ILMath {
47
48        /// <summary>
49        /// Type definitions for possible kernels in kernel ridge regression
50        /// </summary>
51        public enum KRRTypes {
52            /// <summary>
53            /// Linear kernel <c>k(x,y) = x'*y</c>
54            /// </summary>
55            linear,
56            /// <summary>
57            /// Polynominal kernel <c>k(x,y) = (x'*y + c)^d</c>
58            /// </summary>
59            polynomial,
60            /// <summary>
61            /// Gaussian (exponential) kernel <c>k(x,y) = exp(-norm(x-y)/sigma)</c>
62            /// </summary>
63            gaussian
64        }
65
66        internal delegate ILRetArray<double> ILKRRApplyFunc(ILInArray<double> X, ILKRRResult C);
67
68        /// <summary>
69        /// Encapsulates the result of a kernel ridge regression and makes it applicable to new data
70        /// </summary>
71        public class ILKRRResult : ILResult<double> {
72
73            private ILKRRApplyFunc m_applyFunc;
74            internal ILKRRApplyFunc ApplyFunc {
75                get { return m_applyFunc; }
76                set { m_applyFunc = value; }
77            }
78            internal ILArray<double> Shift = returnType<double>();
79            internal ILArray<double> X = returnType<double>();
80            internal KRRTypes Kernel = KRRTypes.gaussian;
81            internal ILArray<double> KernelParameter = returnType<double>();
82            internal ILArray<double> Alpha = returnType<double>();
83
84            public override ILRetArray<double> Apply(ILInArray<double> Data, params ILBaseArray[] arguments)
85            {
86                if (IsDisposed)
87                    throw new InvalidOperationException("this result has already been disposed");
88                return m_applyFunc(Data,this);
89            }
90
91            public override void Dispose() {
92                if (IsDisposed)
93                    return;
94                base.Dispose();
95                if (!object.ReferenceEquals(Shift, null))
96                    Shift.Dispose();
97                if (!object.ReferenceEquals(X, null))
98                    X.Dispose();
99                if (!object.ReferenceEquals(KernelParameter, null))
100                    KernelParameter.Dispose();
101                if (!object.ReferenceEquals(Alpha,null))
102                    Alpha.Dispose();
103            }
104        }
105
106        /// <summary>
107        /// Calculate the kernel ridge regression (KRR) of X to Y
108        /// </summary>
109        /// <param name="X">Input data X, observations in columns</param>
110        /// <param name="Y">Input data Y (regression target)</param>
111        /// <returns>Kernel ridge regression model with linear kernel and no regularization</returns>
112        public static ILResult<double> krr(ILInArray<double> X, ILInArray<double> Y) {
113            return krr(X, Y, KRRTypes.linear, 1, -1);
114        }
115        /// <summary>
116        /// Calculate the kernel ridge regression (KRR) of X to Y
117        /// </summary>
118        /// <param name="X">Input data X, observations in columns</param>
119        /// <param name="Y">Input data Y (regression target)</param>
120        /// <param name="kernel">Kernel type to use</param>
121        /// <param name="Kernelparam">Kernel parameters, scalar constant, if invalid: defaults to 1</param>
122        /// <param name="Regularization">Regularization constant (set to -1 to have no regularization)</param>
123        /// <returns>The kernel ridge regression model</returns>
124        public static ILResult<double> krr(ILInArray<double> X, ILInArray<double> Y, KRRTypes kernel,
125                                           ILInArray<double> Kernelparam, ILInArray<double> Regularization) {
126            using (ILScope.Enter(X, Y, Kernelparam, Regularization)) {
127                ILArray<double> x = check(X);
128                ILArray<double> y = check(Y);
129                ILArray<double> kernelparam = check(Kernelparam, (a) => { return (a.IsScalar) ? a : 1.0; });
130                ILArray<double> regularization = check(Regularization, (a) => { return (a.IsScalar) ? a : 0.0; });
131
132                ILKRRResult C = new ILKRRResult();
133
134                C.ApplyFunc = krrapply;
135                // centering
136                C.Shift.a = mean(x, 1);
137                //% X = X - repmat(C.Shift,1,n);
138                x = x - repmat(C.Shift, 1, x.S[1]);
139                ILArray<double> K = empty<double>();
140                switch (kernel) {
141                    case KRRTypes.linear:
142                        K = multiply(x.T, x);                                                               // linear: K = X'*X;
143                        break;
144                    case KRRTypes.polynomial:
145                        K = pow((multiply(x.T, x) + 1), kernelparam);                                           // % K = (X'*X + 1).^kernelparameter;
146                        break;
147                    case KRRTypes.gaussian:
148                        K = gkernel(x, x, kernelparam);                                                     // % K = gkernel(X,X,kernelparameter);
149                        break;
150                }
151                C.X.a = x.C;
152                C.Kernel = kernel;
153                C.KernelParameter.a = kernelparam;
154               
155                y = (y.S[0] < y.S[1]) ? y.T : y.C;                                                          // if (size(y,1) < size(y,2)) y = y'; end
156                if (regularization == 0) {                                                                    // if (~exist('regularization') || regularization == 0)
157                    ILArray<double> U = empty<double>();                                                      //    [U,V] = eig(K);
158                    ILArray<double> V = empty<double>();
159                    U.a = eigSymm(K, V);
160                    ILArray<double> c = diag(V);                                                              //    c = diag(V); 
161                    ILArray<double> minLooErr = double.PositiveInfinity;                                      //    minLooErr = inf;
162                    //                                                                                        //    cr = find(c > 1e-12)';                                 
163                    foreach (int ci in find(c > 1e-12)) {                                                     //    for ci = cr     
164                        using (ILScope.Enter()) {
165                            ILArray<double> Sd = c / (c + c[ci]);                                             //        Sd = c./(c + c(ci));
166                            ILArray<double> S = multiply(U, diag(Sd), U.T);                                   //        S = U * diag(Sd) * U';
167                            ILArray<double> looer = sum(y - multiply(S, y) / pow(1 - diag(S), 2.0)) / x.S[1];   //        looerr = sum((y - S*y)./(1-diag(S)).^2) ./ n;
168                            if (minLooErr > looer) {                                                          //        if (minLooErr > looerr)
169                                minLooErr.a = looer;                                                          //            minLooErr = looerr;
170                                regularization.a = c[ci];                                                     //            regularization = c(ci);
171                            }                                                                                 //        end
172                        }                                                                                     //    end
173                    }
174                    // since we already have U, lets compute alpha that way
175                    C.Alpha.a = multiply(U, diag(1 / (c + regularization)), U.T, y).T;                        //    C.alpha = U * diag(1./(c + regularization)) * U' * y;
176                    //                                                                                        //    C.alpha = C.alpha';
177                    return C;                                                                                 //    return;
178                } else {                                                                                      //end
179                    C.Alpha.a = linsolve((K + diag(zeros(1, x.S[1]) + regularization)), y).T;                 // C.alpha = (K + diag(zeros(1,n) + regularization))\y;
180                    return C;                                                                                 // C.alpha = C.alpha';
181                }                                                                                             //end
182            }
183        }
184
185        private static ILRetArray<double> gkernel(ILInArray<double> X1, ILInArray<double> X2, ILInArray<double> kernelparam) {
186            using (ILScope.Enter(X1, X2, kernelparam)) {
187                int d = X1.S[0], n1 = X1.S[1], n2 = X2.S[1];                    // [d,n1] = size(X1); n2 = size(X2,2);
188                ILArray<double> XX = multiply(X1.T, X2);                        // XX = X1'*X2;
189                ILArray<double> K = repmat(sum(X2 * X2, 0), n1, 1);             // K = repmat(sum(X2.^2,1),n1,1);
190                K.a = K + repmat(sum(X1 * X1, 0).T, 1, n2);                     // K = K + repmat(sum(X1.^2,1)',1,n2);
191                K.a = K - 2 * XX;                                               // K = K - 2.*XX;
192                K.a = exp(-K / (2 * d * kernelparam * kernelparam));            // K = exp(-K./(2 * d * w.^2));
193                return K;
194            }
195        }
196
197        private static ILRetArray<double> krrapply(ILInArray<double> X, ILKRRResult C) {
198            using (ILScope.Enter(X)) {
199                ILArray<double> x = check(X); 
200                int d = x.S[0], n = x.S[1];                                     // [d,n] = size(X);
201                x.a = x - repmat(C.Shift, 1, n);                                  // X = X - repmat(C.Shift,1,n);
202                ILArray<double> K = empty<double>();
203                switch (C.Kernel) {
204                    case KRRTypes.linear:
205                        K.a = multiply(x.T, C.X);                               // K = X'*C.X;
206                        break;
207                    case KRRTypes.polynomial:
208                        K.a = pow(multiply(x.T, C.X) + 1, C.KernelParameter);    // K = (X'*C.X + 1).^C.parameter;
209                        break;
210                    case KRRTypes.gaussian:
211                        K.a = gkernel(x, C.X, C.KernelParameter);                // K = gkernel(X,C.X,C.parameter);
212                        break;
213                }
214                return sum(C.Alpha * K, 1).T;                              //Y = sum(repmat(C.alpha,n,1).*K,2)'; 
215            }
216        }
217
218    }
219}
Note: See TracBrowser for help on using the repository browser.