# Dependent gaussian processes

I need a help with gaussian processes. I am implementing dependent gaussian processes as on this paper

Boyle, Phillip, and Marcus Frean. "Dependent gaussian processes." Advances in neural information processing systems. 2004.

$P\left( {y,m,\sigma } \right) \sim N(C\left| \sum \right.) % MathType!MTEF!2!1!+- % feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn % hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr % 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq-Jc9 % vqaqpepm0xbba9pwe9Q8fs0-yqaqpepae9pg0FirpepeKkFr0xfr-x % fr-xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaamiuamaabm % aabaGaamyEaiaacYcacaWGTbGaaiilaiabeo8aZbGaayjkaiaawMca % aiablYJi6iaad6eacaGGOaGaam4qamaaeeaabaGaeyyeIuoacaGLhW % oacaGGPaaaaa!44BD!$

with a covariance of the form

${C^y} = {C^u} + \sigma I % MathType!MTEF!2!1!+- % feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn % hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr % 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq-Jc9 % vqaqpepm0xbba9pwe9Q8fs0-yqaqpepae9pg0FirpepeKkFr0xfr-x % fr-xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaam4qamaaCa % aaleqabaGaamyEaaaakiabg2da9iaadoeadaahaaWcbeqaaiaadwha % aaGccqGHRaWkcqaHdpWCcaWGjbaaaa!3E66!$

where C is equal to $${C_u} = \left[ {\begin{array}{*{20}{c}}{{C_{11}}}&{{C_{12}}}\\{{C_{21}}}&{{C_{22}}}\end{array}} \right] % MathType!MTEF!2!1!+- % feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn % hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr % 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq-Jc9 % vqaqpepm0xbba9pwe9Q8fs0-yqaqpepae9pg0FirpepeKkFr0xfr-x % fr-xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaam4qamaaBa % aaleaacaWG1baabeaakiabg2da9maadmaabaqbaeqabiGaaaqaaiaa % doeadaWgaaWcbaGaaGymaiaaigdaaeqaaaGcbaGaam4qamaaBaaale % aacaaIXaGaaGOmaaqabaaakeaacaWGdbWaaSbaaSqaaiaaikdacaaI % XaaabeaaaOqaaiaadoeadaWgaaWcbaGaaGOmaiaaikdaaeqaaaaaaO % Gaay5waiaaw2faaaaa!44CA!$$ and

$$\begin{array}{l}Cov_{ij}^u = \sum\limits_{m = 1}^M {\frac{{{{\left( {2\pi } \right)}^{\frac{D}{2}}}{v_{mi}}{v_{mj}}}}{{\sqrt {\left| {{A_{mj}} + {A_{mi}}} \right|} }}\exp \left( { - \frac{1}{2}{{\left( {{d_s} - \left[ {{u_{mi}} - {u_{mj}}} \right]} \right)}^T}\sum \left( {{d_s} - \left[ {{u_{mi}} - {u_{mj}}} \right]} \right)} \right)} \\\sum = {A_{mi}}{({A_{mi}} + {A_{mj}})^{ - 1}}{A_{mj}}\end{array} % MathType!MTEF!2!1!+- % feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn % hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr % 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq-Jc9 % vqaqpepm0xbba9pwe9Q8fs0-yqaqpepae9pg0FirpepeKkFr0xfr-x % fr-xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGceaqabeaacaWGdb % Gaam4BaiaadAhadaqhaaWcbaGaamyAaiaadQgaaeaacaWG1baaaOGa % eyypa0ZaaabCaeaadaWcaaqaamaabmaabaGaaGOmaiabec8aWbGaay % jkaiaawMcaamaaCaaaleqabaWaaSaaaeaacaWGebaabaGaaGOmaaaa % aaGccaWG2bWaaSbaaSqaaiaad2gacaWGPbaabeaakiaadAhadaWgaa % WcbaGaamyBaiaadQgaaeqaaaGcbaWaaOaaaeaadaabdaqaaiaadgea % daWgaaWcbaGaamyBaiaadQgaaeqaaOGaey4kaSIaamyqamaaBaaale % aacaWGTbGaamyAaaqabaaakiaawEa7caGLiWoaaSqabaaaaOGaciyz % aiaacIhacaGGWbWaaeWaaeaacqGHsisldaWcaaqaaiaaigdaaeaaca % aIYaaaamaabmaabaGaamizamaaBaaaleaacaWGZbaabeaakiabgkHi % TmaadmaabaGaamyDamaaBaaaleaacaWGTbGaamyAaaqabaGccqGHsi % slcaWG1bWaaSbaaSqaaiaad2gacaWGQbaabeaaaOGaay5waiaaw2fa % aaGaayjkaiaawMcaamaaCaaaleqabaGaamivaaaakiabggHiLpaabm % aabaGaamizamaaBaaaleaacaWGZbaabeaakiabgkHiTmaadmaabaGa % amyDamaaBaaaleaacaWGTbGaamyAaaqabaGccqGHsislcaWG1bWaaS % baaSqaaiaad2gacaWGQbaabeaaaOGaay5waiaaw2faaaGaayjkaiaa % wMcaaaGaayjkaiaawMcaaaWcbaGaamyBaiabg2da9iaaigdaaeaaca % WGnbaaniabggHiLdaakeaacqGHris5cqGH9aqpcaWGbbWaaSbaaSqa % aiaad2gacaWGPbaabeaakiaacIcacaWGbbWaaSbaaSqaaiaad2gaca % WGPbaabeaakiabgUcaRiaadgeadaWgaaWcbaGaamyBaiaadQgaaeqa % aOGaaiykamaaCaaaleqabaGaeyOeI0IaaGymaaaakiaadgeadaWgaa % WcbaGaamyBaiaadQgaaeqaaaaaaa!8F4E!$$

and a Likelihood

$L = - \frac{1}{2}\log \left| C \right| - \frac{1}{2}y'{C^{ - 1}}y - \frac{N}{2}\log (2\pi ) % MathType!MTEF!2!1!+- % feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn % hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr % 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq-Jc9 % vqaqpepm0xbba9pwe9Q8fs0-yqaqpepae9pg0FirpepeKkFr0xfr-x % fr-xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaamitaiabg2 % da9iabgkHiTmaalaaabaGaaGymaaqaaiaaikdaaaGaciiBaiaac+ga % caGGNbWaaqWaaeaacaWGdbaacaGLhWUaayjcSdGaeyOeI0YaaSaaae % aacaaIXaaabaGaaGOmaaaacaWG5bGaai4jaiaadoeadaahaaWcbeqa % aiabgkHiTiaaigdaaaGccaWG5bGaeyOeI0YaaSaaaeaacaWGobaaba % GaaGOmaaaaciGGSbGaai4BaiaacEgacaGGOaGaaGOmaiabec8aWjaa % cMcaaaa!51EC!$

My question is with respect to the partial derivative of the likelihood (gradients)in respect to each one of the hyperparameters . Until now, I had been using lsqnonlin in matlab. but fails too many times to obtain a good aproximation to the simple data. I want to use minimize from Rasmussen or other strategies but these optimizations need the gradients. As I need to program for nxm input and outputs how can I code the gradients for each Hyperparameter? any suggestions are welcome.

I am having this results afeter many restarts

I upload my code to GITHUB if anyone want to comment or help me to develop this project are welcome

https://github.com/arizawilmer/MIMO-Gaussian-Process-Identification-Matlab

• Which algorithm options did you use in lsqnonlin? There are probably better algorithm options in some other solvers (not in Optimization Toolbox). You may be able to use something like ADiMat to do the needed gradient, and even potentially Hessian calculations via matrix level automatic differentiation under MATLAB. You can call the gradient (and Hessian if you use Hessian option) in the user function you provide the optimizer to calculate gradient (and Hessian. Use of Hessian in a trust region Newton method will probably be the most robust. Quasi-Newton using gradient might be faster/easier. – Mark L. Stone Oct 5 '16 at 6:34
• For ADiMat, which is free, go to sc.informatik.tu-darmstadt.de/res/adimat/index.en.jsp . – Mark L. Stone Oct 5 '16 at 6:35
• Hi,Thank you I will check and update the question as i get the solution as I just start to code for more than two by two I am using Levenberg-Marquardt I know that is near usable but for speed test of the code was alright until 2X2. There is any form to share the code for comments of the community? – Wilmer Ariza Oct 5 '16 at 22:27
• I upload my code to github if anyone want to make suggestion, feel free to do it. github.com/arizawilmer/… – Wilmer Ariza Oct 6 '16 at 0:39
• I also read the paper. Do you find the source code for the paper? The only question that puzzles me is that I can't realize the 3D figure in the paper. Could you give me a hand? – yang Apr 9 '18 at 11:44

The derivatives can be calculated as $$\frac{{\partial \mathcal{L}\left( \Theta \right)}}{{\partial {\Theta _i}}} = - \frac{1}{2}trace\left( {{K^{ - 1}}\frac{{\partial K}}{{\partial {\Theta _i}}}} \right) + \frac{1}{2}{y^T}{K^{ - 1}}\frac{{\partial K}}{{\partial {\Theta _i}}}{K^{ - 1}}y$$
$$\frac{{\partial K}}{{\partial {\Theta _i}}}$$ can be calculated depending on the position on the matrix there are 4 different positions for a 4x4 matrix that gives us two kernels $${k_{i = j}}$$ and $${k_{i \ne j}}$$ each kernel is defined by 3 and 4 hyperparameters