Proper Complex Gaussian Processes for Regression
Rafael~Boloix-Tortosa*
F.~Javier~Pay\'an-Somet
Eva Arias-de-Reyna
and Juan~Jos\'e~Murillo-Fuentes% <-this % stops a space \thanks{R. Boloix-Tortosa
F.~Javier~Pay\'an-Somet
Eva Arias-de-Reyna
and Juan~Jos\'e~Murillo-Fuentes are with the Department of Signal Theory and Communications
University of Seville
Spain. e-mail: rboloix@us.es. Phone: +34 954488132.
....lastname \thanks{...
Complex-valued signals are used in the modeling of many systems in engineering and science, hence being of fundamental interest. Often, random complex-valued signals are considered to be proper. A proper complex random variable or process is uncorrelated with its complex conjugate. This assumption is a good model of the underlying physics in many problems, and simplifies the computations. While linear processing and neural networks have been widely studied for these signals, the development of complex-valued nonlinear kernel approaches remains an open problem. In this paper we propose Gaussian processes for regression as a framework to develop 1) a solution for proper complex-valued kernel regression and 2) the design of the reproducing kernel for complex-valued inputs, using the convolutional approach for cross-covariances. In this design we pay attention to preserve, in the complex domain, the measure of similarity between near inputs. The hyperparameters of the kernel are learned maximizing the marginal likelihood using Wirtinger derivatives. Besides, the approach is connected to the multiple output learning scenario. In the experiments included, we first solve a proper complex Gaussian process where the cross-covariance does not cancel, a challenging scenario when dealing with proper complex signals. Then we successfully use these novel results to solve some problems previously proposed in the literature as benchmarks, reporting a remarkable improvement in the estimation error. % % while Also, the assumption of proper complex random variables and processes being uncorrelated with their complex conjugate simplifies the % %They as they allow to process a pair of signals Lately the nonlinear kernel based regression is paper we investigate Gaussian processes for regression (GPR) for complex-valued signals. We develop a novel complete statistical description of GPR predictions for proper complex-valued signals, along with the optimization procedure to set hyperparameters, applying Wirtinger derivatives. %Hence, the formulation is constrained to provide proper outputs. %Involved matrix operations benefit from this new formulation, since complex matrix operations, such as matrix inversion, are faster and more accurate than real counterparts in double dimension. Also, as a by-product, how to yield complex covariance functions from real ones is also explained from this framework, while preserving the intuitive physical meaning of a measure of similarity between near inputs. Finally, the Wirtinger derivatives are used to find the hyperparameters that maximize the marginal likelihood. % % \cite{Ehrlich70, \cite{Golub96
arxiv2015.tex

%% bare_jrnl.tex
%% V1.4
%% 2012/12/27
%% by Michael Shell
%% see http://www.michaelshell.org/
%% for current contact information.
%%
%% This is a skeleton file demonstrating the use of IEEEtran.cls
%% (requires IEEEtran.cls version 1.8 or later) with an IEEE journal paper.
%%
%% Support sites:
%% http://www.michaelshell.org/tex/ieeetran/
%% http://www.ctan.org/tex-archive/macros/latex/contrib/IEEEtran/
%% and
%% http://www.ieee.org/



% *** Authors should verify (and, if needed, correct) their LaTeX system  ***
% *** with the testflow diagnostic prior to trusting their LaTeX platform ***
% *** with production work. IEEE's font choices can trigger bugs that do  ***
% *** not appear when using other class files.                            ***
% The testflow support page is at:
% http://www.michaelshell.org/tex/testflow/


%%*************************************************************************
%% Legal Notice:
%% This code is offered as-is without any warranty either expressed or
%% implied; without even the implied warranty of MERCHANTABILITY or
%% FITNESS FOR A PARTICULAR PURPOSE! 
%% User assumes all risk.
%% In no event shall IEEE or any contributor to this code be liable for
%% any damages or losses, including, but not limited to, incidental,
%% consequential, or any other damages, resulting from the use or misuse
%% of any information contained here.
%%
%% All comments are the opinions of their respective authors and are not
%% necessarily endorsed by the IEEE.
%%
%% This work is distributed under the LaTeX Project Public License (LPPL)
%% ( http://www.latex-project.org/ ) version 1.3, and may be freely used,
%% distributed and modified. A copy of the LPPL, version 1.3, is included
%% in the base LaTeX documentation of all distributions of LaTeX released
%% 2003/12/01 or later.
%% Retain all contribution notices and credits.
%% ** Modified files should be clearly indicated as such, including  **
%% ** renaming them and changing author support contact information. **
%%
%% File list of work: IEEEtran.cls, IEEEtran_HOWTO.pdf, bare_adv.tex,
%%                    bare_conf.tex, bare_jrnl.tex, bare_jrnl_compsoc.tex,
%%                    bare_jrnl_transmag.tex
%%*************************************************************************

% Note that the a4paper option is mainly intended so that authors in
% countries using A4 can easily print to A4 and see how their papers will
% look in print - the typesetting of the document will not typically be
% affected with changes in paper size (but the bottom and side margins will).
% Use the testflow package mentioned above to verify correct handling of
% both paper sizes by the user's LaTeX system.
%
% Also note that the "draftcls" or "draftclsnofoot", not "draft", option
% should be used if it is desired that the figures are to be displayed in
% draft mode.
%
%\documentclass[journal]{IEEEtran}
\documentclass[11pt,onecolumn] {IEEEtran} 
%\documentclass[10pt,twocolumn,twoside] {IEEEtran}
%
% If IEEEtran.cls has not been installed into the LaTeX system files,
% manually specify the path to it like:
% \documentclass[journal]{../sty/IEEEtran}

% Some very useful LaTeX packages include:
% (uncomment the ones you want to load)


% *** MISC UTILITY PACKAGES ***
%
%\usepackage{ifpdf}
% Heiko Oberdiek's ifpdf.sty is very useful if you need conditional
% compilation based on whether the output is pdf or dvi.
% usage:
% \ifpdf
%   % pdf code
% \else
%   % dvi code
% \fi
% The latest version of ifpdf.sty can be obtained from:
% http://www.ctan.org/tex-archive/macros/latex/contrib/oberdiek/
% Also, note that IEEEtran.cls V1.7 and later provides a builtin
% \ifCLASSINFOpdf conditional that works the same way.
% When switching from latex to pdflatex and vice-versa, the compiler may
% have to be run twice to clear warning/error messages.


% *** CITATION PACKAGES ***
%
%\usepackage{cite}
% cite.sty was written by Donald Arseneau
% V1.6 and later of IEEEtran pre-defines the format of the cite.sty package
% \cite{} output to follow that of IEEE. Loading the cite package will
% result in citation numbers being automatically sorted and properly
% "compressed/ranged". e.g., [1], [9], [2], [7], [5], [6] without using
% cite.sty will become [1], [2], [5]--[7], [9] using cite.sty. cite.sty's
% \cite will automatically add leading space, if needed. Use cite.sty's
% noadjust option (cite.sty V3.8 and later) if you want to turn this off
% such as if a citation ever needs to be enclosed in parenthesis.
% cite.sty is already installed on most LaTeX systems. Be sure and use
% version 4.0 (2003-05-27) and later if using hyperref.sty. cite.sty does
% not currently provide for hyperlinked citations.
% The latest version can be obtained at:
% http://www.ctan.org/tex-archive/macros/latex/contrib/cite/
% The documentation is contained in the cite.sty file itself.


% *** GRAPHICS RELATED PACKAGES ***
%
\ifCLASSINFOpdf
   \usepackage[pdftex]{graphicx}
  % declare the path(s) where your graphic files are
  % \graphicspath{{../pdf/}{../jpeg/}}
  % and their extensions so you won't have to specify these with
  % every instance of \includegraphics
  % \DeclareGraphicsExtensions{.pdf,.jpeg,.png}
\else
  % or other class option (dvipsone, dvipdf, if not using dvips). graphicx
  % will default to the driver specified in the system graphics.cfg if no
  % driver is specified.
  % \usepackage[dvips]{graphicx}
  % declare the path(s) where your graphic files are
  % \graphicspath{{../eps/}}
  % and their extensions so you won't have to specify these with
  % every instance of \includegraphics
  % \DeclareGraphicsExtensions{.eps}
\fi
% graphicx was written by David Carlisle and Sebastian Rahtz. It is
% required if you want graphics, photos, etc. graphicx.sty is already
% installed on most LaTeX systems. The latest version and documentation
% can be obtained at: 
% http://www.ctan.org/tex-archive/macros/latex/required/graphics/
% Another good source of documentation is "Using Imported Graphics in
% LaTeX2e" by Keith Reckdahl which can be found at:
% http://www.ctan.org/tex-archive/info/epslatex/
%
% latex, and pdflatex in dvi mode, support graphics in encapsulated
% postscript (.eps) format. pdflatex in pdf mode supports graphics
% in .pdf, .jpeg, .png and .mps (metapost) formats. Users should ensure
% that all non-photo figures use a vector format (.eps, .pdf, .mps) and
% not a bitmapped formats (.jpeg, .png). IEEE frowns on bitmapped formats
% which can result in "jaggedy"/blurry rendering of lines and letters as
% well as large increases in file sizes.
%
% You can find documentation about the pdfTeX application at:
% http://www.tug.org/applications/pdftex


% *** MATH PACKAGES ***
%
\usepackage[cmex10]{amsmath}
\interdisplaylinepenalty=2500
\usepackage{amssymb}
%:Paquete para notaciones específicas
\usepackage{notacion2}
\newtheorem{definition}{Definition}[section]

\ifx\pdfoutput\undefined
% we are running LaTeX, not pdflatex
\usepackage{graphicx}
\else
% we are running pdflatex, so convert .eps files to .pdf
\usepackage[pdftex]{graphicx}
%\usepackage{graphicx}
\usepackage{epstopdf}
%\usepackage{epspdfconversion}
\fi
%%%%%%%%%%%%%%%%%%%%%%%
% Bloque para que las figuras se generen en el texto
\usepackage[svgnames, x11names]{xcolor}
%\usepackage{dtsc-creafig}

% Definición de las dimensiones de la figura
%\newlength\figureheight 
%\newlength\figurewidth 

\usepackage{environ}


\NewEnviron{equationw}{
\begin{equation}\begin{split}
\BODY
\end{split}\end{equation}
}



% A popular package from the American Mathematical Society that provides
% many useful and powerful commands for dealing with mathematics. If using
% it, be sure to load this package with the cmex10 option to ensure that
% only type 1 fonts will utilized at all point sizes. Without this option,
% it is possible that some math symbols, particularly those within
% footnotes, will be rendered in bitmap form which will result in a
% document that can not be IEEE Xplore compliant!
%
% Also, note that the amsmath package sets \interdisplaylinepenalty to 10000
% thus preventing page breaks from occurring within multiline equations. Use:
%\interdisplaylinepenalty=2500
% after loading amsmath to restore such page breaks as IEEEtran.cls normally
% does. amsmath.sty is already installed on most LaTeX systems. The latest
% version and documentation can be obtained at:
% http://www.ctan.org/tex-archive/macros/latex/required/amslatex/math/


% *** SPECIALIZED LIST PACKAGES ***
%
%\usepackage{algorithmic}
% algorithmic.sty was written by Peter Williams and Rogerio Brito.
% This package provides an algorithmic environment fo describing algorithms.
% You can use the algorithmic environment in-text or within a figure
% environment to provide for a floating algorithm. Do NOT use the algorithm
% floating environment provided by algorithm.sty (by the same authors) or
% algorithm2e.sty (by Christophe Fiorio) as IEEE does not use dedicated
% algorithm float types and packages that provide these will not provide
% correct IEEE style captions. The latest version and documentation of
% algorithmic.sty can be obtained at:
% http://www.ctan.org/tex-archive/macros/latex/contrib/algorithms/
% There is also a support site at:
% http://algorithms.berlios.de/index.html
% Also of interest may be the (relatively newer and more customizable)
% algorithmicx.sty package by Szasz Janos:
% http://www.ctan.org/tex-archive/macros/latex/contrib/algorithmicx/


% *** ALIGNMENT PACKAGES ***
%
%\usepackage{array}
% Frank Mittelbach's and David Carlisle's array.sty patches and improves
% the standard LaTeX2e array and tabular environments to provide better
% appearance and additional user controls. As the default LaTeX2e table
% generation code is lacking to the point of almost being broken with
% respect to the quality of the end results, all users are strongly
% advised to use an enhanced (at the very least that provided by array.sty)
% set of table tools. array.sty is already installed on most systems. The
% latest version and documentation can be obtained at:
% http://www.ctan.org/tex-archive/macros/latex/required/tools/


% IEEEtran contains the IEEEeqnarray family of commands that can be used to
% generate multiline equations as well as matrices, tables, etc., of high
% quality.


% *** SUBFIGURE PACKAGES ***
%\ifCLASSOPTIONcompsoc
  \usepackage[caption=false,font=normalsize,labelfont=sf,textfont=sf]{subfig}
%\else
%  \usepackage[caption=false,font=footnotesize]{subfig}
%\fi
% subfig.sty, written by Steven Douglas Cochran, is the modern replacement
% for subfigure.sty, the latter of which is no longer maintained and is
% incompatible with some LaTeX packages including fixltx2e. However,
% subfig.sty requires and automatically loads Axel Sommerfeldt's caption.sty
% which will override IEEEtran.cls' handling of captions and this will result
% in non-IEEE style figure/table captions. To prevent this problem, be sure
% and invoke subfig.sty's "caption=false" package option (available since
% subfig.sty version 1.3, 2005/06/28) as this is will preserve IEEEtran.cls
% handling of captions.
% Note that the Computer Society format requires a larger sans serif font
% than the serif footnote size font used in traditional IEEE formatting
% and thus the need to invoke different subfig.sty package options depending
% on whether compsoc mode has been enabled.
%
% The latest version and documentation of subfig.sty can be obtained at:
% http://www.ctan.org/tex-archive/macros/latex/contrib/subfig/


% *** FLOAT PACKAGES ***
%
%\usepackage{fixltx2e}
% fixltx2e, the successor to the earlier fix2col.sty, was written by
% Frank Mittelbach and David Carlisle. This package corrects a few problems
% in the LaTeX2e kernel, the most notable of which is that in current
% LaTeX2e releases, the ordering of single and double column floats is not
% guaranteed to be preserved. Thus, an unpatched LaTeX2e can allow a
% single column figure to be placed prior to an earlier double column
% figure. The latest version and documentation can be found at:
% http://www.ctan.org/tex-archive/macros/latex/base/


%\usepackage{stfloats}
% stfloats.sty was written by Sigitas Tolusis. This package gives LaTeX2e
% the ability to do double column floats at the bottom of the page as well
% as the top. (e.g., "\begin{figure*}[!b]" is not normally possible in
% LaTeX2e). It also provides a command:
%\fnbelowfloat
% to enable the placement of footnotes below bottom floats (the standard
% LaTeX2e kernel puts them above bottom floats). This is an invasive package
% which rewrites many portions of the LaTeX2e float routines. It may not work
% with other packages that modify the LaTeX2e float routines. The latest
% version and documentation can be obtained at:
% http://www.ctan.org/tex-archive/macros/latex/contrib/sttools/
% Do not use the stfloats baselinefloat ability as IEEE does not allow
% \baselineskip to stretch. Authors submitting work to the IEEE should note
% that IEEE rarely uses double column equations and that authors should try
% to avoid such use. Do not be tempted to use the cuted.sty or midfloat.sty
% packages (also by Sigitas Tolusis) as IEEE does not format its papers in
% such ways.
% Do not attempt to use stfloats with fixltx2e as they are incompatible.
% Instead, use Morten Hogholm'a dblfloatfix which combines the features
% of both fixltx2e and stfloats:
%
% \usepackage{dblfloatfix}
% The latest version can be found at:
% http://www.ctan.org/tex-archive/macros/latex/contrib/dblfloatfix/



%\ifCLASSOPTIONcaptionsoff
%  \usepackage[nomarkers]{endfloat}
% \let\MYoriglatexcaption\caption
% \renewcommand{\caption}[2][\relax]{\MYoriglatexcaption[#2]{#2}}
%\fi
% endfloat.sty was written by James Darrell McCauley, Jeff Goldberg and 
% Axel Sommerfeldt. This package may be useful when used in conjunction with 
% IEEEtran.cls'  captionsoff option. Some IEEE journals/societies require that
% submissions have lists of figures/tables at the end of the paper and that
% figures/tables without any captions are placed on a page by themselves at
% the end of the document. If needed, the draftcls IEEEtran class option or
% \CLASSINPUTbaselinestretch interface can be used to increase the line
% spacing as well. Be sure and use the nomarkers option of endfloat to
% prevent endfloat from "marking" where the figures would have been placed
% in the text. The two hack lines of code above are a slight modification of
% that suggested by in the endfloat docs (section 8.4.1) to ensure that
% the full captions always appear in the list of figures/tables - even if
% the user used the short optional argument of \caption[]{}.
% IEEE papers do not typically make use of \caption[]'s optional argument,
% so this should not be an issue. A similar trick can be used to disable
% captions of packages such as subfig.sty that lack options to turn off
% the subcaptions:
% For subfig.sty:
% \let\MYorigsubfloat\subfloat
% \renewcommand{\subfloat}[2][\relax]{\MYorigsubfloat[]{#2}}
% However, the above trick will not work if both optional arguments of
% the \subfloat command are used. Furthermore, there needs to be a
% description of each subfigure *somewhere* and endfloat does not add
% subfigure captions to its list of figures. Thus, the best approach is to
% avoid the use of subfigure captions (many IEEE journals avoid them anyway)
% and instead reference/explain all the subfigures within the main caption.
% The latest version of endfloat.sty and its documentation can obtained at:
% http://www.ctan.org/tex-archive/macros/latex/contrib/endfloat/
%
% The IEEEtran \ifCLASSOPTIONcaptionsoff conditional can also be used
% later in the document, say, to conditionally put the References on a 
% page by themselves.



% *** PDF, URL AND HYPERLINK PACKAGES ***
%
%\usepackage{url}
% url.sty was written by Donald Arseneau. It provides better support for
% handling and breaking URLs. url.sty is already installed on most LaTeX
% systems. The latest version and documentation can be obtained at:
% http://www.ctan.org/tex-archive/macros/latex/contrib/url/
% Basically, \url{my_url_here}.




% *** Do not adjust lengths that control margins, column widths, etc. ***
% *** Do not use packages that alter fonts (such as pslatex).         ***
% There should be no need to do such things with IEEEtran.cls V1.6 and later.
% (Unless specifically asked to do so by the journal or conference you plan
% to submit to, of course. )


% correct bad hyphenation here



\begin{document}
%
% paper title
% can use linebreaks \\ within to get better formatting as desired
% Do not put math or special symbols in the title.
\title{Proper Complex Gaussian Processes for Regression}
%
%
% author names and IEEE memberships
% note positions of commas and nonbreaking spaces ( ~ ) LaTeX will not break
% a structure at a ~ so this keeps an author's name from being broken across
% two lines.
% use \thanks{} to gain access to the first footnote area
% a separate \thanks must be used for each paragraph as LaTeX2e's \thanks
% was not built to handle multiple paragraphs
%

\author{Rafael~Boloix-Tortosa*,
        F.~Javier~Pay\'an-Somet,
        Eva Arias-de-Reyna,
        and Juan~Jos\'e~Murillo-Fuentes% <-this % stops a space
\thanks{R. Boloix-Tortosa, F.~Javier~Pay\'an-Somet, Eva Arias-de-Reyna,
        and Juan~Jos\'e~Murillo-Fuentes are with the Department of Signal Theory and Communications, University of Seville, Spain. e-mail: rboloix@us.es. Phone: +34 954488132.}% <-this % stops a space
%\thanks{J. Doe and J. Doe are with Anonymous University.}% <-this % stops a space
%\thanks{Thanks to Spanish government (Ministerio de Educaci\'on y Ciencia, TEC2012-38800-C03-02) and European Union (FEDER) for funding.}% <-this % stops a space
}

% note the % following the last \IEEEmembership and also \thanks - 
% these prevent an unwanted space from occurring between the last author name
% and the end of the author line. i.e., if you had this:
% 
% \author{....lastname \thanks{...} \thanks{...} }
%                     ^------------^------------^----Do not want these spaces!
%
% a space would be appended to the last name and could cause every name on that
% line to be shifted left slightly. This is one of those "LaTeX things". For
% instance, "\textbf{A} \textbf{B}" will typeset as "A B" not "AB". To get
% "AB" then you have to do: "\textbf{A}\textbf{B}"
% \thanks is no different in this regard, so shield the last } of each \thanks
% that ends a line with a % and do not let a space in before the next \thanks.
% Spaces after \IEEEmembership other than the last one are OK (and needed) as
% you are supposed to have spaces between the names. For what it is worth,
% this is a minor point as most people would not even notice if the said evil
% space somehow managed to creep in.



% The paper headers
%\markboth{IEEE Signal Processing Letters,~Vol.~X, No.~X, December~201X}%
%{Shell \MakeLowercase{\textit{et al.}}: Complex Proper Gaussian Processes for Regression}
% The only time the second header will appear is for the odd numbered pages
% after the title page when using the twoside option.
% 
% *** Note that you probably will NOT want to include the author's ***
% *** name in the headers of peer review papers.                   ***
% You can use \ifCLASSOPTIONpeerreview for conditional compilation here if
% you desire.




% If you want to put a publisher's ID mark on the page you can do it like
% this:
%\IEEEpubid{0000--0000/00\$00.00~\copyright~2012 IEEE}
% Remember, if you use this you must call \IEEEpubidadjcol in the second
% column for its text to clear the IEEEpubid mark.



% use for special paper notices
%\IEEEspecialpapernotice{(Invited Paper)}


\hyphenation{op-tical net-works semi-conduc-tor hy-per-pa-ra-me-ters ge-ne-ra-tion ge-ne-ra-ted}

% make the title area
\maketitle

% As a general rule, do not put math, special symbols or citations
% in the abstract or keywords.
\begin{abstract}
Complex-valued signals are used in the modeling of many systems in engineering and science, hence being of fundamental interest. Often, random complex-valued signals are considered to be proper. A proper complex random variable or process is uncorrelated with its complex conjugate. This assumption is a good model of the underlying physics in many problems, and simplifies the computations. While linear processing and neural networks have been widely studied for these signals, the development of complex-valued nonlinear kernel approaches remains an open problem. In this paper we propose Gaussian processes for regression as a framework to develop 1) a solution for proper complex-valued kernel regression and 2) the design of the reproducing kernel for complex-valued inputs, using the convolutional approach for cross-covariances. In this design we pay attention to preserve, in the complex domain, the measure of similarity between near inputs. The hyperparameters of the kernel are learned maximizing the marginal likelihood using Wirtinger derivatives. Besides, the approach is connected to the multiple output learning scenario. In the experiments included, we first solve a proper complex Gaussian process where the cross-covariance does not cancel, a challenging scenario when dealing with proper complex signals. Then we successfully use these novel results to solve some problems previously proposed in the literature as benchmarks, reporting a remarkable improvement in the estimation error.


%
% while   Also, the assumption of proper complex random variables and processes being uncorrelated with their complex conjugate simplifies the 
%
%They as they allow to process a pair of signals Lately  the nonlinear kernel based regression is paper we investigate Gaussian processes for regression (GPR) for complex-valued signals. We develop a novel complete statistical description of GPR predictions for proper complex-valued signals, along with the optimization procedure to set hyperparameters, applying Wirtinger derivatives. %Hence, the formulation is constrained to provide proper outputs. 
%Involved matrix operations benefit from this new formulation, since complex matrix operations, such as matrix inversion, are faster and more accurate than real counterparts in double dimension. Also, as a by-product, how to yield complex covariance functions from real ones is also explained from this framework, while preserving the intuitive physical meaning of a measure of similarity between near inputs. Finally, the Wirtinger derivatives are used to find the hyperparameters that maximize the marginal likelihood. 
%
% \cite{Ehrlich70}, \cite{Golub96}

 \end{abstract}

% Note that keywords are not normally used for peerreview papers.
\begin{IEEEkeywords}
Gaussian processes, regression, proper complex processes, kernel methods.
\end{IEEEkeywords}






% For peer review papers, you can put extra information on the cover
% page as needed:
% \ifCLASSOPTIONpeerreview
% \begin{center} \bfseries EDICS Category: 3-BBND \end{center}
% \fi
%
% For peerreview papers, this IEEEtran command inserts a page break and
% creates the second title. It will be ignored for other modes.
\IEEEpeerreviewmaketitle

\setlength{\arraycolsep}{0.1em} % para que no deje mucha separación en las ecuaciones

\section{Introduction}
% The very first letter is a 2 line initial drop letter followed
% by the rest of the first word in caps.
% 
% form to use if the first word consists of a single letter:
% \IEEEPARstart{A}{demo} file is ....
% 
% form to use if you need the single drop letter followed by
% normal text (unknown if ever used by IEEE):
% \IEEEPARstart{A}{}demo file is ....
% 
% Some journals put the first two words in caps:
% \IEEEPARstart{T}{his demo} file is ....
% 
% Here we have the typical use of a "T" for an initial drop letter
% and "HIS" in caps to complete the first word.


\IEEEPARstart{C}{omplex-valued} signals model a vast range of nowadays systems in science and engineering such as telecommunications, optics, electromagnetics, and acoustics among others. The main advantage of using complex-valued signals is the availability of processing the real and imaginary parts as a single signal. For these reasons, complex-valued signal processing is of fundamental interest.

Signal processing for complex-valued signals has been widely studied in the linear case, see \cite{Schreier06} and references therein. A proper complex random signal is uncorrelated with its complex conjugate \cite{Neeser93}. Properness is useful because it simplifies computations and in many cases models the underlying physics of the problem at hand. Most of the linear solutions in the real case admit a proper complex signals version with minor changes \cite{Tulay11}.  For the processing of improper signals further development of the linear solutions are needed \cite{Tulay11,Schreier06}.

The nonlinear processing of complex-valued signals has been addressed from the point of view of neural networks, \cite{hirose13} and, recently, using reproducing kernel Hilbert spaces (RKHS) \cite{Scholkopf02}. Some complex kernels have been lately proposed for classification  \cite{Steinwart06}, regression \cite{OgunfunmiP11,Bouboulis12,Tobar12} and mainly for kernel principal component analysis \cite{Papaioannou14}. %and real 
%The kernel based solutions proposed for regression adopts real case structures. 
Regarding regression, in \cite{Bouboulis11} the authors propose a complex-valued kernel based in the results in \cite{Steinwart06} and face the derivative of cost functions by using Wirtinger's derivatives. Same kernel is adopted in \cite{OgunfunmiP11}. As discussed later in this paper the resulting approach involves properness. Besides, the kernel used is neither stationary nor isotropic, and it may suffer from numerical problems. In \cite{Tobar12} the authors review the kernel design to improve the previous solution with a kernel they denote as independent. The resulting kernel yields also proper complex-valued outputs. The kernel is stationary, but again it is not isotropic in the complex-valued input space, as the real and imaginary parts of the input are split and fed to different real valued kernels. Hence, these previous works do not report results for isotropic and stationary kernels that may better model the underlying physics of some systems. Also, the structure of the kernel remains more rigid than needed. These drawbacks make these solutions not powerful enough to learn a wide range of systems. 
%
The previous approaches have been developed in the framework of kernel least squares and for adaptive filtering, although the extension to batch processing is immediate. We bring here Gaussian processes for regression (GPR) as a useful tool to better study this problem of kernel complex-valued nonlinear regression \cite{Boloix14}. 
%In \cite{Boloix14} we gave a first step by rewriting the GPR model to solve a complex-valued linear system in digital communications. Here, we first develop in  \SEC{GPRC} a full statistical formulation of the complex GPR %(GPRC) 
%
%We outlined some of  this approach in \cite{} to solve the inversion of linear mixtures of digital communication signals,   

GPRs are well known techniques in machine learning \cite{Rasmussen06,PerezCruz13gp}. They have been successfully applied to regression, classification and dimensionality reduction. GPs can be interpreted as a family of kernel methods with the additional advantage of providing a full conditional statistical description for the predicted variable, where hyperparameters can be learned by maximizing the marginal likelihood, avoiding cross-validation. 
%
Properties and features of random complex-valued Gaussian processes have been widely studied \cite{Miller69}, facilitating the extension of GPR for the complex-valued case proposed in this paper.  The benefits are twofold. On the one hand we have a novel tool for nonlinear regression providing probabilistic outputs. Second, since this Bayesian approach is connected to the regularization point of view \cite{Alvarez12}, we prove the structure of the reproducing kernel and study the conditions to be valid for proper complex outputs. 
%
The solution of GPR in the general complex-valued case is quite involved, since the pseudo-covariance is to be taken into account. In the proper scenario the pseudo-covariance cancels, and the formulation of Gaussian processes for real-valued signals can be adapted \cite{Neeser93}. We describe this new solution, providing the mean and the covariance of the complex-valued outputs. 
%
We particularly focus on the kernel, to improve previous solutions. We derive its structure, being complex valued, where the real part of the kernel is given by the covariance of the real part of the outputs plus the covariance of the imaginary part of the outputs, while the imaginary part of the kernel describes the cross-covariance between real and imaginary parts of the outputs. We prove that the real and imaginary parts of the kernel can be designed with different features. But we conclude that the imaginary part, in addition to be skew-symmetric, must be constructed to ensure the whole covariance to be semi-definite positive, i.e. a reproducing kernel or covariance matrix \cite{Taubock12}. 
%
We also pay attention to the modeling of physical systems, and propose a kernel for the real part based on the exponential kernel. At this point, when measuring similarity between inputs we resort to the absolute value of the complex difference between inputs. As a result, the kernel is isotropic and stationary. Other real kernels may be adapted using this metric, to better model the physics of the problem at hand. 
The construction of the imaginary part is challenging. Conditions on this matrix to ensure properness and to be a valid kernel quite limit the degree of freedom in the kernel design, and the systems we may model. We resort to the convolution approach \cite{Boyle05,Calder07} to ensure the whole matrix to be a reproducing kernel. We propose a model for the cross-covariance that explains a positive and negative correlation of the real part of the output with the imaginary one for a positive and negative delay, canceling at origin.   
%
%
Every system modeled with complex-valued signals can be rewritten with real-valued ones by splitting the output signal into a composite two-dimensional vector, with real and imaginary parts. For the sake of comparison, the multiple output learning (MOL) \cite{Boyle05,Alvarez12} is proposed to solve this composite kernel complex-valued regression problem and adapt it to the proper complex case. In this scenario we lose the intuition about complex values, and we cannot deal with inputs and outputs easily pairing real and imaginary parts to measure similarity between them. Also, we may lose some of the computational benefits of working with complex numbers,  
such as in the computation of the matrix inversion,  \cite{Ehrlich70}, \cite{Golub96}. %which are faster and more accurate than real counterparts in double dimension \cite{Ehrlich70}, \cite{Golub96}. 
However, we prove that both approaches yield equivalent results, gaining some further intuition on the kernel properties and structure.   

One of the advantages of GPR is the availability of estimating hyperparameters by maximizing the marginal likelihood. Here we face the derivation of a function of complex-valued parameters and variables. We cope with this problem by using Wirtinger's derivatives. 
%
To illustrate the good performance of this approach, we include a first experiment where we learn a Gaussian process with unknown hyperparameters, where we assume the cross-covariance does not cancel. Then we compare to the algorithms in \cite{Bouboulis12,Tobar12} using the same scenario than in the experiments described in \cite{Bouboulis12} as benchmark. Due to the improvements derived from the proposals in this paper, we achieve remarkable good results. At this point we conjecture that whenever we do not know the cause or model of the cross-covariance or this cross-covariance is negligible, assuming a real-valued reproducing kernel is a good option.

%A Gaussian process is a collection of random variables, any finite number of which have a joint real Gaussian distribution. 

%Recently, a few kernel methods have been proposed for complex-valued nonlinear regression \cite{Bouboulis11, Bouboulis12, Tobar12}. These works mainly focus on the kernel design and the least mean square algorithm. However, to the best of our knowledge, no full description of Gaussian processes (GPs) regression for complex-valued signals has been developed so far. 
%GPR for real-valued signals is well known. 

%In this paper we investigate the complex GPR (CGPR) formulation for proper complex signals. In \cite{Boloix14} we gave a first step by rewriting the GPR model to solve a complex-valued linear system in digital communications. Here, we first develop in  \SEC{GPRC} a full statistical formulation of the complex GPR %(GPRC) 
%predictions for proper complex-valued observations.

The rest of the paper is organized as follows. Next section includes the derivation of GPR for proper complex-valued signals. In \SEC{Cov} we exploit this formulation to study the structure and construction of the kernel, compare it to the kernel in the MOL approach, and provide some examples for the kernel. %with the resulting from the MOL its structure and define the procedure to derive complex kernels from well known real ones. Finally, 
We develop in \SEC{ML} the optimization procedure to set the kernel hyperparameters applying Wirtinger's calculus and patterned complex-valued matrix derivatives \cite{Hjorungnes08}. The good performance of this approach is illustrated in \SEC{Exp}, learning a unknown process and comparing to previous approaches in the channel equalization problem. We end with conclusions and an appendix.
%where proper CGPR is used to learn a given proper complex-valued Gaussian process with unknown hyperparameters. And we compare to previously We conclude remarking that the regression is efficient, in the sense that the formulation takes into account the properness of the  output. 
%Another benefit of working with this novel formulation is that it allows complex matrix operations, such as matrix inversion, which are faster and more accurate than real counterparts in double dimension \cite{Ehrlich70}, \cite{Golub96}. 
%To summarize, the main contributions of this letter are: a) the development of a complete statistical description of GPR predictions for proper complex-valued signal; b) the generic extension of real kernels or covariance functions to complex ones, allowing well known real valued kernels to be adapted to the complex case preserving the intuitive physical meaning of a measure of similarity for near inputs; c) model hyperparameters optimization using Wirtinger’s Calculus and patterned complex-valued matrix derivatives.


The notations used in the paper are as follows. If $\vect{A}$ is a matrix, $\entry{\vect{A}}{l}{q}$ is its $(l,q)$ entry. For a vector, $\vect{a}$, $a_l$ denotes its $l$-th entry.
%To denote a range of values from $r$ to $q$, $r\!:\!q$ is used, e.g., $\entry{\vect{A}}{r:q}{k}$ is a column vector, the values $r$ to $q$ of the column $k$ of $\vect{A}$. $\entry{\vect{A}}{:}{k}$ is the $k$th column of $\vect{A}$. 
$\vect{A}\trs$ represents the transpose of $\vect{A}$ and $\vect{A}\her$ its Hermitian. % and $\vect{A}\pinv$ its psudoinverse. 
$\trace{\vect{A}}$ is the trace of $\vect{A}$, and $\det{\vect{A}}$ its determinant. % and $\diag{\vect{A}}$ is a column vector with the diagonal of $\vect{A}$.  
To denote the $i$-th sample we use $(i)$, both for a vector, $\vect{a}(i)$, and a scalar, $a(i)$. %The $r$-th entry of vector ${\vect{a}}$ is referred to as ${\vect{a}}(r)$. 
${\vect{a}}_n$ is a vector of $n$ entries while $\vect{A}_n$ is a matrix with $n$ columns. The real and imaginary parts are denoted by subindex $\textrm{r}$ and $\textrm{j}$, respectively, i.e. $\vect{a}\rr=\Re(\vect{a})$ and $\vect{a}\jj=\Im(\vect{a})$, and throughout the text $\textrm{j}=\sqrt{-1}$. To denote the complex Gaussian distribution with mean vector $\boldsymbol{\mu}$, covariance matrix $\K$ and pseudo-covariance matrix $\matr{\tilde{K}}$ we use $\calg{N}\left(\boldsymbol{\mu},\K,\matr{\tilde{K}}\right)$.
%$\diag{\vect{a}}$ is a diagonal matrix whose diagonal is vector $\vect{a}$. 
%$||\vect{a}||$ is the $L_2$ norm of vector $\vect{a}$.  %denotes an $\nsc$ module cyclic shift 
%$\vect{\vect{A}}$ is a column vector with the columns of $\vect{A}$, from left to right, stacked downwards.  
$\E[{\cdot}]$ denotes statistical expectation. %Finally, $\hdp$ and $\hdd$ stand for the Hadamard (pointwise) product and division, respectively.
%$\vect{A}$. 
%and $\vect{A}\pinv$ is its pseudoinverse.
% and the entries of $\vect{A}\cnj$ are the complex conjugate of the elements of 




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Proper Complex Gaussian Process Regression}\LABSEC{GPRC}
%GPs are characterized by their mean and covariance functions, denoted by $\mu(\x)$ and $\k(\x_{i},\x_{l})$
%The output of the regression, i.e. the latent function in the model in \EQ{regression}, can be evaluated 
GP for regression can be presented as a nonlinear regressor that expresses the input-output relation through function $f(\x)$, known as latent function, that follows a GP and underlies the regression problem
\begin{equation}\LABEQ{regression}
{\y}=\f(\x)+\epsilon,
\end{equation}
where the error, $\epsilon$, in the estimation of the output, $\y$, is modeled as additive zero-mean Gaussian noise.
GPs are the natural nonlinear Bayesian extension to the linear minimum mean-squared error (LMMSE) and Wiener filtering \cite{MurilloGP09,PerezCruz13gp}, thus, GPs provide the correct approach to solve an MMSE filter nonlinearly. 

For any input set $\{ (\x(i)) | i = 1, . . . , \n\}$, where $\x({i}) \in {\CN}^d$ are complex-valued column vectors of dimension $\d$, the latent function in \EQ{regression} can be designed and evaluated to provide a multidimensional
Gaussian complex-valued random vector $\fv_\n=\left[{f}(\x({1})), . . . , {f}(\x(\n))\right]^{\top}$, where $f(\x(i))\in \CN$. In this paper we focus on both complex-valued input and output. The simpler real-valued input and complex-valued output case can be easily solved from the results herein. 
The immediate approach when dealing with complex-valued signals is to consider the real and imaginary components as separate real signals \cite{PerezCruz08b}, for both inputs and outputs. If we wish real and imaginary parts to be considered dependent, we may resort to multiple output Gaussian processes to relate them \cite{Boyle05,Alvarez12}. However, in these approaches we lose the complex value notion,  specially for the input.

%Hence, real GPR in double input dimension can be applied. 

A complex Gaussian vector is characterized not only by its mean vector $\boldsymbol{\mu}$ and covariance matrix $\K= \mathbb{E}\left[(\fv_\n-\boldsymbol{\mu})(\fv_\n-\boldsymbol{\mu})\her\right]$, but also by its complementary covariance or pseudo-covariance matrix $\matr{\tilde{K}}=\mathbb{E}\left[(\fv_\n-\boldsymbol{\mu})(\fv_\n-\boldsymbol{\mu})^\top\right]$, \cite{Schreier06}. We use the notation $\calg{N}\left(\boldsymbol{\mu},\K,\matr{\tilde{K}}\right)$ for a joint Gaussian probability distribution to emphasize this fact.
%:
%\begin{align}
%p(\fv_\n|\X_\n)&=\calg{N}\left(\boldsymbol{\mu},\K,\matr{\tilde{K}}\right)\nonumber\\
%&=\frac{1}{\pi^n \det^{1/2}\aug{\K}}\exp{\left(-\frac{1}{2}(\aug{\fv}_\n-\aug{\boldsymbol{\mu}})\her\aug{\K}\inv(\aug{\fv}_\n-\aug{\boldsymbol{\mu}})\right)},
%\end{align}
%where $\X_\n=\left[\x({1}), . . . , \x(\n)\right]$, $\aug{\fv}_\n=\left[{\fv}^\top_\n \;{\fv}\her_\n\right]^\top$ is the augmented vector for $\fv_\n$, $\aug{\boldsymbol{\mu}}=\left[\boldsymbol{\mu}^\top \;\boldsymbol{\mu}\her\right]^\top$ is the augmented mean vector of $\fv_\n$, and $\aug{\K}$ is the augmented covariance matrix of $\fv_\n$ given by
%\begin{equation}
%\aug{\K}=\left[\begin{array}{c c}
%\K & \psc{\K}\\
%\psc{\K}^* & \K^*\\
%\end{array}\right].
%\end{equation}
%Therefore, in order to characterize a complex GP, the mean function, $\mu(\x)$, and both the covariance and complementary covariance functions, $\k(\x(i),\x(l))=(\K)_{il}$ and $\psc{\k}(\x(i),\x(l))=(\psc{\K})_{il}$ respectively, are needed. %Here $(\K)_{il}$ is the $(i,l)$ entry of $\K$. 
When the complementary covariance matrix is zero, a complex Gaussian random vector is regarded as \emph{proper} \cite{Neeser93,Ollila08}. %Also, a proper complex Gaussian process has zero complementary covariance function. 
Without loss of generality, we consider here a zero-mean process for the latent function \cite{Neeser93}. By imposing properness, $(\psc{\K})_{il}=\psc{\k}(\x(i),\x(l))=0$ for all $i,l$. The prior for this proper complex Gaussian process yields %$p(\fv_\n|\X_\n)=\calg{N}\left(\vect{0},\K,\matr{0}\right)$.
\begin{equationw}\label{eq:propergaussian}
p(\fv_\n|\X_\n)=\calg{N}\left(\vect{0},\K,\matr{0}\right)=\frac{1}{\pi^n \det{\K}}\exp{\left(-{\fv}_\n\her{\K}\inv{\fv}_\n\right)},
\end{equationw}
where $\X_\n=\left[\x({1}), . . . , \x(\n)\right]$, and $({\K})_{il}={\k}(\x(i),\x(l))$ for all  $i,l$ in the input set, being ${\k}(\x(i),\x(l))$ the covariance function.
% and the last zero in the parameters of the normal distribution denotes a null pseudo-covariance matrix.

In the learning process we condition the output of the GPR for some new observations %, ${\X\new }_{m}=\left[\x\new({1}), . . . , \x\new(m)\right]$, 
given the training set $\tset=\{(\x(i),y({i})) | i = 1, . . . , \n\}=\{\X_\n,\yv_\n\}$, where the outputs  $\yv_\n=\left[{\y}({1}), . . . , {\y}(\n)\right]^{\top}$, $\y(i)\in \CN$, for a given set of observations $\X_\n$, are known. We can compute this posterior through the joint distribution of the to-be-estimated outputs and the training data $\yv_\n$. This distribution can be constructed by using \EQ{propergaussian} and the marginal likelihood, $p\left(\yv_\n| \X_\n \right)$, as follows. %as the product of the prior and the likelihood divided by the marginal likelihood (or evidence). 
%Bayes' rule gives the posterior density for $\fv_\n$,
%\begin{equation}
%p\left(\fv_\n|\mathcal{D}_\n\right)=\frac{p\left(\yv_\n|\fv_\n \right)p\left(\fv_\n|\X_\n \right)}{p\left(\yv_\n|\X_\n \right)},
%\end{equation}
%where $p\left(\yv_\n|\fv_\n \right)$ is the likelihood for the latent function at the training set and $p\left(\yv_\n|\X_\n \right)$ is the marginal likelihood or evidence of the model. 
We assume that the additive noise $\epsilon$ in \EQ{regression} follows an i.i.d. proper (circular) complex Gaussian distribution with zero mean and variance $\sigma^2_\epsilon$.
The samples in the training set are i.i.d., hence the likelihood for the latent function at the training set, $p\left(\yv_\n|\fv_\n \right)$, is given by a factorized model
\begin{equation}
p\left(\yv_\n| \fv_\n \right)= \prod_{i=1}^n {p(y({i})|f(\x({i})))},
\end{equation}
where $p(y({i})|f(\x({i})))=\calg{N}(f(\x({i})),\sigma^2_\epsilon,0)$.
%, i.e., 
%\begin{equation}
%p(y_i)=\frac{1}{\pi \det^{1/2}\aug{\matr{\Sigma}}_{\epsilon}}e^{\left\{-\frac{1}{2}(\aug{y}_i-\aug{f}(\x_i))\her\aug{\matr{\Sigma}}_{\epsilon}\inv(\aug{y}_i-\aug{f}(\x_i))\right\}},
%\end{equation}
%where $\aug{y}_i=\left[{y}_i \;{y}^*_i\right]^\top$, $\aug{f}(\x_i)=\left[{f}(\x_i)) \;{f}(\x_i))^*\right]^\top$, and
%\begin{equation}
%\aug{\matr{\Sigma}}_{\epsilon}=\left[\begin{array}{c c}
%{\Sigma}_{\epsilon} & {\tilde\Sigma}_{\epsilon}\\
%{\tilde\Sigma}_{\epsilon}^* & {\Sigma}_{\epsilon}\\
%\end{array}\right]=\sigma^2_\epsilon\left[\begin{array}{c c}
%1 & \rho\\
%\rho^* & 1\\
%\end{array}\right].
%\end{equation}
Therefore, the likelihood is a proper complex multidimensional Gaussian $p\left(\yv_\n| \fv_\n \right)= \calg{N}\left(\fv_\n,\sigma^2_\epsilon\id{\n},\matr{0}\right)$. This likelihood and the prior in \EQ{propergaussian} yield the marginal likelihood or evidence 
\begin{equation}\label{eq:marlikelihood}
p\left(\yv_\n| \X_\n \right)=\int{p\left(\yv_\n| \fv_\n \right)p(\fv_\n|\X_\n)d\fv_\n} =\calg{N}\left(\vect{0},\C,\matr{0}\right),
\end{equation}
where  $\C=\K+\sigma_\epsilon^2\id{\n}$.
Note that $\yv_\n$ is also proper Gaussian. Furthermore, $\yv_\n$ and $\fv_\n$ are cross-proper, as the complementary cross-covariance matrix $\mathbb{E}\left[\yv_\n\fv_\n^\top\right]=\matr{0}$.
Hence, $\yv_\n$ and $\fv_\n$ are jointly proper \cite{Schreier06}, i.e., the composite complex random vector $\left[\yv_\n^\top,\fv_\n^\top\right]^\top$ is Gaussian proper. If we now have a number $m$ of test input vectors ${\X\new }_{m}=\left[\x\new({1}), . . . , \x\new({m})\right]$, the joint distribution of the training outputs $\yv_\n$ and ${\fv\new}_{m}=\left[{f}(\x\new({1})), . . . , {f}(\x\new({m}))\right]^{\top}$ is as follows: 
\begin{equation}\label{eq:yf}
\left[ \begin{array}{c}
\yv_\n\\
{\fv\new}_{m}\\
\end{array}\right]\sim\calg{N}\left( \matr{0}, \left[ \begin{array}{ccc}
\C & {} &\K(\X_\n,{\X\new }_{m})\\
\K({\X\new }_{m},\X_\n) & {} & \K({\X\new }_{m},{\X\new }_{m})\\
\end{array}\right],\vect{0}\right),
\end{equation}
where the entry $(i,l)$ of matrix $\K(\X_\n,{\X\new }_{m})$ is $\k(\x({i}),\x\new({l}))$, and in a similar way are defined $\K({\X\new }_{m},\X_\n)$ and $\K({\X\new }_{m},{\X\new }_{m})$. From these results, the estimated probabilistic output is the conditional distribution of ${\fv\new}_{m}$ given $\yv_\n$:
\begin{equation} \label{eq:posteriorf}
p\left({\fv\new}_{m}| {\X\new }_{m},\X_\n,\yv_\n\right)=\calg{N}\left(\boldsymbol{\mu}_{{\fv\new}_{m}},\matr{\Sigma}_{{\fv\new}_{m}},\matr{0}\right),
\end{equation}
where 
\begin{align}
\boldsymbol{\mu}_{{\fv\new}_{m}}&=\K({\X\new }_{m},\X_\n)\C\inv\yv_\n, \LABEQ{mediaCPGP}\\
\matr{\Sigma}_{{\fv\new}_{m}}&=\K({\X\new }_{m},{\X\new }_{m})-\K({\X\new }_{m},\X_\n)\C\inv\K(\X_\n,{\X\new }_{m}).\label{eq:varCPGP}
\end{align}
%In the case that there is only one test input vector $\x\new$, 
%\begin{align}
%{\mu}_{\f\new} &= \kv\new\her\C\inv\yv_\n,\label{eq:mediaCPGPsingle}\\
%{\Sigma}_{\f\new}&={k}(\x\new,\x\new)-\kv\new\her\C\inv\kv\new,
%\end{align}
%with $\kv\new=\left[\k(\x_1,\x\new),\k(\x_2,\x\new),...,\k(\x_\n,\x\new)\right]^\top$.


The predictive distribution is also proper Gaussian. %, with a mean given by the posterior mean of the weights multiplied by the test input. 
The predictive covariance is a quadratic form of the test and training inputs, showing that the predictive uncertainties grow with the magnitude of the covariance matrices involved, as one would expect for a linear model.  
%Note that if the covariance between test and training input data vanishes, the term to the right in \EQ{varCPGP} vanishes and the covariance is that of the prior in \EQ{propergaussian}.



%\begin{tikzpicture}[auto, node distance=2.5cm,>=latex,line width=1.2pt, color=black]
%	\node [input] (entrada) {};
%	\node [buffer gate US,draw, inner sep=2pt,minimum height=1cm,right of=entrada] (gain) {$a$};
%	\node[circle, draw, inner sep=0pt,minimum size=1cm,right of=gain](sumador) {};
%	\draw[] (node cs:name=sumador,anchor=north) -- (node cs:name=sumador,anchor=south);
%      \draw[] (node cs:name=sumador,anchor=west) -- (node cs:name=sumador,anchor=east);
%      %
%	\node[bloque, right of=sumador] (sistema3) {Detector};	                                             
%	\node [output, right of=sistema3] (output) {};
%	\node [input, below of=sumador] (ruido) {};
%	%
%	\draw [->] (entrada) -- node[name=x] {$s_{i}$} (gain);
%	\draw[->] (gain) -- node[name=x2] {}  (sumador);
%	\draw [->] (ruido) -- node[name=N] {$N$} (sumador);
%	\draw [->] (sumador) -- node[name=r] {$r$} (sistema3);
%	\draw [->] (sistema3) -- node [name=y] {$\hat{s}$}(output);	
%\end{tikzpicture}

%\vspace*{-2.3cm}%\hspace{8cm}%
%\begin{tikzpicture}[ scale=1.8]%, show background rectangle]
%\draw (2,2) circle (1);
%\draw (1 mm, 10 pt) -- (4 em, 1);

%\path[clip] (1,1)--(2,2)--(3,1)--cycle;
%\path[fill=blue!50] (2, 1.7) circle (.8);
%\path[draw,fill=green!50] (1, 1.7) circle (.8);
%
%\draw (0, 0) rectangle (2, 1);
%\draw[color=red] (0, 0) circle (.5);






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Complex Covariance Functions}\LABSEC{Cov}
%
The covariance function is a key tool in GPR: it encodes our assumptions about the function that we wish to learn and measures similarity between inputs. There are some well known examples of covariance functions used for real-valued applications, such as the squared exponential or Gaussian covariance function and the inhomogeneous polynomial kernel, among others. These kernels assume that training points that are near to a test point should be informative about the prediction at that point \cite{Rasmussen06}. The definition of a kernel for complex-valued GPR should provide the same measure of similarity between the inputs. 

%We propose a general method to make use of well known real-valued examples of covariance functions in complex GPR. 
In the proper complex-valued Gaussian processes framework we can derive the kernel following a more principled approach as follows.
Consider a zero-mean complex Gaussian vector $\fv_\n={\fv_\n}\rr+\j{\fv_\n}\jj$, with ${\fv_\n}\rr$ its real part and ${\fv_\n}\jj$ its imaginary part. The covariance matrix $\K= \mathbb{E}\left[\fv_\n\fv_\n\her\right]$ is \cite{Schreier06}: 
\begin{equation}\label{eq:compleskernel}
\K=\K\rrrr+\K\jjjj+\j\left(\K\jjrr-\K\rrjj\right),
\end{equation}
where $\K\rrrr$ and $\K\jjjj$ $\in\RN_+^{n\times n}$ are the covariance matrices of real and imaginary parts of $\fv_\n$, respectively, and $\K\rrjj= \mathbb{E}\left[{\fv_\n}\rr{\fv_\n}\jj^\top\right]=\K\jjrr\trs$ $\in\RN^{n\times n}$ is the cross-covariance matrix of real and imaginary parts. At this point it is most important to notice that $\K$ must be a valid covariance matrix for the Gaussian process. Real-valued covariance functions can now be used to write out the three real covariance matrices $\K\rrrr$, $\K\jjjj$ and $\K\rrjj$, but they are interrelated. Given two inputs, the complex covariance function of the process is composed by three real covariance functions, 
\begin{equation}\label{eq:complexCV}
\k(\x,\x')=k\rrrr(\x,\x')+k\jjjj(\x,\x')+\j\left(k\rrjj(\x',\x)-k\rrjj(\x,\x')\right).
\end{equation}
%which take into account the covariances and cross-covariance of real and imaginary parts of the process. 


For the particular case of a proper complex Gaussian vector, $\K\rrrr=\K\jjjj$ and $\K\jjrr=\K\rrjj\trs=-\K\rrjj$ \cite{Schreier06}. Therefore, in the case of a proper complex Gaussian process, functions $k\rrrr(\x,\x')=k\jjjj(\x,\x')$, and $k\rrjj(\x,\x')$ must yield either a null or skew-symmetric cross-covariance matrix $\K\rrjj$. The covariance matrix yields
\begin{equation}\LABEQ{complesProperkernel}
\K=\K\rr +\j\K\jj,
\end{equation}
where $\K\rr =2\K\rrrr$ and $\K\jj=-2\K\rrjj$. Following the guidelines in \cite{Alvarez12} we may conclude that the proposed kernel in \EQ{complesProperkernel} is a valid reproducing kernel if it is a covariance matrix for the Bayesian Gaussian process framework. Hence, it must be a Hermitian positive semi-definite matrix \cite{Taubock12}. $\K\rrrr$ must be a symmetric and positive semi-definite matrix, since the marginals of the joint probability functions of the real and imaginary parts must be also covariance matrices as discussed later in this section. It follows that the condition $\vect{v}\her \K \vect{v}\geq0$ for any $\vect{v}\in\CN$ yields,
\begin{align}\LABEQ{kercon}
\vect{v}\her \K \vect{v}=\vect{v}\rr\trs \K\rr \vect{v}\rr+\vect{v}\jj\trs \K\rr \vect{v}\jj-2 \vect{v}\rr\trs \K\jj \vect{v}\jj \geq 0
\end{align}
where the first terms to the right of the equality are greater or equal to cero, since $\K\rr$ is positive semi-definite and we used $\vect{v}\rr\trs \K\jj \vect{v}\rr=$$\vect{v}\jj\trs \K\jj \vect{v}\jj=0$ since $\K\jj$ is skew-symmetric. 
%and $j\vect{v}\jj\trs \K\rr \vect{v}\jj=j$\vect{v}\jj\trs \K\rr \vect{v}\rr
In the following subsection we analyze this problem by decomposing it into the real and imaginary parts.

%%%
%n=4;
%B=randn(n,n); Tr=triu(B); B=Tr-Tr';
%A=randn(n,n); A=A*A';
%C=[A B;B' A]; eig(C)
%v=randn(n,1)+sqrt(-1)*randn(n,1); w=[real(v); imag(v)];
%[v'*(A+sqrt(-1)*B')*v w'*C*w]
%%ojo a los por dos porque las pdfs se definen de otra forma
%[v'*((A+sqrt(-1)*B')*2)^-1*v w'*(C^-1)*w/2]
%real(v)'*A*real(v)+imag(v)'*A*imag(v)-real(v)'*B'*imag(v)+imag(v)'*B'*real(v)
%Kr=A;Ki=B';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Complex GPR as multiple output GPR}\LABSEC{MO}

%Following the guidelines in \cite{Alvarez12}, the proposed reproducing kernel in \EQ{complesProperkernel} could be proved to be valid if is semidefinite precisely because it is defined within the Bayesian framework of Gaussian processes.
%
% In this subsection we This can also be proved by In this subsection gain a better insight in this we relate it to the real valued multiple output case to provide gain a better insight in the properties of this kernel. 

The complex case may be tackled by mapping the complex value into a two dimensional vector with real and imaginary parts, as in \cite{PerezCruz08}. Then two GPs can be learned, one for the real part and another for the imaginary part of the output. Either independently or using a MOL or vector scheme \cite{Micchelli05,Boyle05,Alvarez12}. The definition of RKHS valued functions in MOL parallels the one in the scalar, where now the reproducing kernel is matrix valued. In our problem and for the training set, we have scalar reproducing kernels $\K\rrrr$, $\K\rrjj$, $\K\jjrr$ and $\K\jjjj$, that are block matrices in the vector kernel matrix
\begin{align} \label{eq:Cmo}
\K_{\Rext} = 
\left[\begin{array}{c c}
\K\rrrr(\X_\n,\X_\n) & \K\rrjj(\X_\n,\X_\n) \\
\K\jjrr(\X_\n,\X_\n) & \K\rrrr(\X_\n,\X_\n)
\end{array}\right].
\end{align} 

For one single test input, we write the mean of the multiple outputs as follows
\begin{align} \label{eq:mimo}
\left[\begin{array}{c}
\hat{\y}\rr \\ \hat{\y}\jj
\end{array}\right]
= 
\left[\begin{array}{c c}
\kv\rrrr({\x\new},\X_\n) & \kv\rrjj({\x\new},\X_\n) \\
\kv\jjrr({\x\new},\X_\n) & \kv\jjjj({\x\new},\X_\n) 
\end{array}\right]
%*
%\left[\begin{array}{c c}
%\C\rrrr & \C\rrjj \\
%\C\jjrr & \C\jjjj
%\end{array}\right]\inv
\C_{\Rext}\inv
%*
\left[\begin{array}{c}
{\yv}\rr \\ {\yv}\jj
\end{array}\right]
\end{align}
where $\C_{\Rext}=\K_{\Rext}+\sigma_\epsilon^2\id{2\n}$ is the covariance matrix in the multiple output or vector formulation.

If we assume the output to be proper the model above yields
\begin{align} \label{eq:mimopr}
\left[\begin{array}{c}
\hat{\y}\rr \\ \hat{\y}\jj
\end{array}\right]
= 
\left[\begin{array}{c c}
\kv\rrrr({\x\new},\X) & \kv\rrjj({\x\new},\X) \\
-\kv\rrjj({\x\new},\X) & \kv\rrrr({\x\new},\X) 
\end{array}\right]
%*
\left[\begin{array}{c c}
\C\rrrr & \C\rrjj \\
-\C\rrjj & \C\rrrr
\end{array}\right]\inv
%*
\left[\begin{array}{c}
{\yv}\rr \\ {\yv}\jj
\end{array}\right].
\end{align}
where, given that the vector covariance matrix must be symmetric, $-\C\rrjj=\C\rrjj^T$. This formulation yields the solution in \EQ{mediaCPGP} with the kernel as proposed in \EQ{complesProperkernel}. Similar conclusions can be drawn for the covariance of the output in \EQ{varCPGP}. This result for the multiple output case links the resulting kernel in the complex value domain to those defined in the MOL case. Since matrix $\C_{\Rext}$
%\begin{align} \label{eq:Cmo}
%\C_{\Rext} = 
%\left[\begin{array}{c c}
%\C\rrrr & \C\rrjj \\
%-\C\rrjj & \C\rrrr
%\end{array}\right]
%\end{align}
is the covariance matrix of the multiple output, it must be Hermitian and positive semi-definite, i.e., equation \EQ{kercon} must hold true for any $\vect{v}$. Also, since the marginal for the real or the imaginary part of the output is also a Gaussian process, matrices $\K\rrrr$ and $\K\rrjj$ must also be symmetric and positive semi-definite. %It also  proving from the multiple output point of view the conditions for \EQ{complesProperkernel} to be a kernel in the proper case. 

If, in addition,
real and imaginary parts are uncorrelated for any pair of inputs, $\C\rrjj=0$, and \EQ{mimopr} yields
\begin{align} \LABEQ{propind}
\left[\begin{array}{c}
\hat{\y}\rr \\ \hat{\y}\jj
\end{array}\right]
= 
%\left[\begin{array}{c c}
%\kv\rrrr({\X\new},\X)\C\rrrr\inv & 0 \\
%0 & \kv\jjjj({\X\new},\X) \C\jjjj\inv
%\end{array}\right]
%%*
\left[\begin{array}{c}
\kv\rrrr({\x\new},\X)\C\rrrr\inv{\yv}\rr \\ \kv\rrrr({\x\new},\X)\C\rrrr\inv{\yv}\jj
\end{array}\right].
\end{align}  
This result can be easily derived from \EQ{mediaCPGP} with $\K=\K\rr=2\K\rrrr$ in \EQ{complesProperkernel}. Here, \EQ{kercon} always holds true if $\K\rrrr$ is positive semi-definite.%, except for the definition of the kernel and covariance matrices. 

The multiple output formulation is equivalent to the complex-valued one. However, the notion of complex signals is lost. Also, if we split the inputs into real and imaginary parts, the design of the covariance matrix becomes a harder task if we want to measure similarity between complex signals.  In contrast, the complex formulation for GPR presented above allows for a kernel design in a natural way, better managing the complex nature of the output and the inputs, as follows.


%One possible solution is to stack the real and imaginary parts of every input into a $2\d$ vector. But we may easily lose the notion of complex number to deal with the imaginary part as any other dimension. 
%
%, withoubidimensional the dimensions of the input in a different way that the imaginary part. 
%A drawback of this multiple output formulation is that it hinders the development of suitable kernels, when dealing with the complex inputs. %First we have to decide if managing the imaginary part as a new input or a new dimension. Where the last case 
%If in the real case we have a $\d\times\n$ matrix $\X$ with training inputs, we may either think of wither a $\d\times 2\n$ or a $2\d\times \n$  matrix in the multiple output case, where the case $2\d\times \n$ allows to write a solution similar to the proper complex GPR. 
%We should not transform the input 
 
  

%, the training set is of size $2n$ where the $n$ imaginary values of the inputs are stacked after the real values,
%%\begin{align}
%$\X_{2\n}=\left[{\x\rr({1})}, . . . , \x\rr(\n), {\x\jj({1})}, . . . , \x\jj(\n)\right]$
%%\end{align}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Design of Covariance Matrices}
%\subsubsection{Null cross-covariance}


In \cite{Steinwart06,Bouboulis11} it is proposed a complex-valued Gaussian kernel as an extension of the real Gaussian kernel: 
\begin{align}\LABEQ{RK}
&\k_{\CN}(\x,\x')=\exp\left(-({\x-\x'^{*})}^\top(\x-\x'^{*})/\gamma\right)\nonumber\\
%&=\exp\left(-(|\x\rr-\x\rr'|^2-|\x\jj+\x\jj'|^2)/\gamma\right)\nonumber\\
%&\cdot\exp\left(-2\j(\x\rr-\x\rr')\trs(\x\jj+\x\jj'))/\gamma\right)\nonumber\\
%&=\exp(-(|\x\rr-\x\rr'|^2-|\x\jj+\x\jj'|^2)/\gamma)\nonumber\\
&=\exp\left(-(|\x\rr-\x\rr'|^2/\gamma\right)\exp\left(|\x\jj+\x\jj'|^2/\gamma\right)\cdot\left(\cos(2(\x\rr-\x\rr')\trs(\x\jj+\x\jj')/\gamma)-\j\sin(2(\x\rr-\x\rr')\trs(\x\jj+\x\jj')/\gamma)\right),
%=&(\x\rr-\x\rr')\trs(\x\rr-\x\rr')-(\x\jj+\x\jj')\trs(\x\jj+\x\jj')\nonumber\\&+2\j(\x\rr-\x\rr')\trs
\end{align}
where $\x=\x\rr+\j\x\jj$, $\x'=\x'\rr+\j\x'\jj$.
%However, this kernel does not have the intuitive physical meaning of a measure of ``similarity'' of the samples. This fact is remarked 
In \cite{Tobar12}, the authors propose the so-called independent kernel to improve the previous one: 
\begin{equation}\LABEQ{indepCGK}
\k_{ind}(\x,\x')=\kappa_{\Rext}\left(\x\rr,\x'\rr\right)+\kappa_{\Rext}\left(\x\jj,\x'\jj\right)+\j\left(\kappa_{\Rext}\left(\x\rr,\x'\jj\right)-\kappa_{\Rext}\left(\x\jj,\x'\rr\right)\right),
\end{equation}
where $\kappa_{\Rext}$ is a real kernel of real inputs. 

%Note that the previously proposed kernels \EQ{RK} in \cite{Bouboulis11} and \EQ{indepCGK} in \cite{Tobar12} and (\ref{eq:complexCV}) exhibit a different structure. In \EQ{RK} we have a unique matrix with no imaginary part. Furthermore, the similarity between complex-valued inputs is 
%These kernels exhibit a different structure to the proposed one, in \EQ{complexCV}. 
In \EQ{RK} we have that the value of the kernel for $(\x,\x')$ is the complex conjugate of the kernel for $(\x',\x)$. This corresponds to assuming the output is proper with non-null skew-symmetric cross-covariance matrix. %The exponent of the kernel is not a metric, since it is not symmetric: the value of the kernel for $(\x,\x')$ is the complex conjugate of the kernel for $(\x',\x)$. It can be easily proved that it assumes that the output is proper with non-null skew-symmetric cross-covariance matrix. %, and that the same exponential kernel is used for the real and imaginary part of the kernel. 
This kernel measures similarities between real parts while measures dissimilarity between imaginary ones and it is not stationary. It also has an oscillatory behavior. In addition, the exponent in the kernel may easily grow large and positive, causing numerical problems, as later discussed in the experiments. These characteristics are not very useful when modeling the underlying physics of many systems.

The kernel in \EQ{indepCGK} solves the measure of similarity between inputs by checking for the real and imaginary parts independently. At this point, we lose the intuition about the complex nature of inputs. The kernel assumes again that the output process to model is proper complex valued, where the imaginary part is non-null and skew-symmetric. One of the main drawbacks of this kernel is that it is not isotropic, due to the way real and imaginary parts of the inputs have been split in the kernel. For example, if a real exponential kernel is used in \EQ{indepCGK} as proposed in \cite{Tobar12}, $\kappa_{\Rext}=\alpha\exp(-|\x-\x'|/\beta)$ for some hyperparameters $\alpha$ and $\beta$, whenever two inputs are distant enough the kernel vanishes except for similar imaginary parts. For any $\x\jj=\x\jj'$, the covariance yields the maximum value for $\kappa_{\Rext}$, $\k_{ind}(\x,\x')=\alpha$. In the imaginary part of the kernel we have a similar behavior. %The imaginary part is skew-symmetric, but the real and imaginary parts of the output at a given input point are correlated with the imaginary and real parts of the output, respectively, along the whole real and imaginary axis in the input space, with origin at that point, except for the input points around the origin. 
%These properties of the kernel prevents modeling a wide range of physical systems. 

%, the real part of the kernel, measuring the covariance of the real and imaginary parts, 

% Again, this may not be the case in many systems. 

%a unique matrix, that due to the is usually complex. with no imaginary part. Hence cross-covariance is assumed to cancel and the covariance of the real and imaginary parts of the outputs to be equal. However, its main drawback is that the measurement of similarity between inputs, that plays a fundamental role in a Hilbert space, is not properly addressed. Note that the argument 
The proposed kernel in \EQ{complexCV} much better adapts to the problem at hand. We may design the real and the imaginary parts of the kernel with different structure. And if the cross-covariance is null or negligible we may set the imaginary part of the kernel to zero. 
% If we have different behavior for the covariances of the real, and imaginary, parts from the cross-covariance, we may use different kernels functions for $\K\rr$ and $\K\jj$. 
The conditions for these kernel functions to form a reproducing kernel are given in \EQ{kercon}. On the other hand, we propose as measurement of similarity or metric the inner product, $\x\her\x$, used in the complex value space as follows.

\subsubsection{Real kernels} 

Assume the cross-covariance cancels and we resort to \EQ{propind}. We do use a \emph{proper} complex Gaussian process to model the system. But we later condition to the training samples. The resulting output within the input range defined by the training inputs may behave, locally, quite differently from a proper Gaussian process with null cross-covariance. Note that in \EQ{propind} we have the same matrix multiplying the real and imaginary parts of the outputs. If these parts are dependent, the GPR will translate this dependence to the output. Furthermore, if the real part of the output has a different variance than the imaginary part, \EQ{propind} properly scales the estimated output, by using the variance of the training outputs. For these reasons, we conjecture that if we know little about the cross-covariance between real and imaginary parts or it is negligible, using \EQ{propind} is a good choice as long as the complex nature of the inputs is properly addressed.

As already discussed, we propose using the inner product of the inputs, $\x\her\x$, a simple metric in complex numbers, to cope with isotropy. The Gaussian kernel yields
\begin{equation} \LABEQ{expkernel}
\k_{G}(\x,\x')=\exp\left(-(\x-\x')\her(\x-\x')/\gamma\right),
\end{equation}
that it is stationary. Therefore, a suitable choice for %any of the three real covariance functions, 
$k\rrrr(\x,\x')=k\jjjj(\x,\x')$ 
%$k\rrrr(\x,\x')$, $k\jjjj(\x,\x')$ or $k\rrjj(\x',\x)$, 
can be $\k_{G}(\x,\x')$ in \EQ{expkernel} while $k\rrjj(\x,\x')$ could be set to zero if we assume that the cross-covariance between real and imaginary parts cancels. If any further information on the physics of the model is known, other kernels \cite{Rasmussen06} can be derived using this metric. 
%We will use this function in the experiments section.  
%

\subsubsection{Complex-valued kernels}The design of the imaginary part of the kernel in \EQ{complexCV} is a challenging problem that remains open, as it is quite system dependent. Here we propose just one kernel, but the procedure followed could help in other designs. The design should met condition \EQ{kercon}. This condition is not constructive, i.e. it does not help designing the kernel. But it encodes intuitive facts such as the maximum absolute value of the cross-covariance being lower or equal to the maximum absolute value of the covariance. On the other hand, the kernel must be able to explain the dependencies between real and imaginary parts of the output, if known. And we are restricted to skew-symmetric matrices. 
%If we know a dependence exists between real and imaginary parts of the output, we should include it in the kernel design. 
%At this point, it is interesting to note the following. 

We propose to model a system where the real and imaginary parts of the outputs are correlated for delayed points with delay $\mu \in \CN^\d$. Since the covariance matrix must be skew symmetric, the correlation will be positive (or negative) for a delay of $\mu$ and negative (or positive) for a delay of $-\mu$. We must ensure that the kernel corresponds to a covariance matrix. We bring here the convolution approach \cite{Boyle05,Calder07}. We model the output of the process as the filtering of two independent white Gaussian noises, $S\rr$ and $S\jj$, and compute the kernels from the filter responses. This way we met the condition of being a covariance matrix. The filters are designed to model the proposed system, fulfilling the imaginary part of the kernel being skew-symmetric. The convolution process is sketched in Fig. 1 where we use exponential responses for the filters,
\begin{align}\LABEQ{fil}
h\rr(\x)&=h\jj(\x)=v\rr\exp(-(\x)\her(\x)/\gamma)\nonumber\\
h\rrjj(\x)&=v\rrjj\exp(-(\x-\mu)\her(\x-\mu)/\gamma)\nonumber\\
h\jjrr(\x)&=-v\rrjj\exp(-(\x+\mu)\her(\x+\mu)/\gamma)
\end{align}
and $W$ is some zero-mean stationary Gaussian noise with variance $\sigma_W^2$. For these filters, it can be proved, see the Appendix, that up to a multiplying constant the kernel in \EQ{complexCV} yields,
\begin{align}\LABEQ{kerfil}
\k(\x, \x')&=(v\rr^{2}+v\rrjj^{2})\exp\left(-\frac{\mathbf{d}_\x\her\mathbf{d}_\x}{2\gamma}\right)\nonumber\\
&+\j v\rr v\rrjj\left(\exp\left(-\frac{(\mathbf{d}_\x-\mu)\her(\mathbf{d}_\x-\mu)}{2\gamma}\right)\exp\left(-\frac{(\mathbf{d}_\x+\mu)\her(\mathbf{d}_\x+\mu)}{2\gamma}\right)\right).
\end{align}
where $\mathbf{d}_\x=\x'-\x$. Note that the constants $v\rrjj$ $\in \RN$, $v\rr$ $\in \RN$ and $\mu$ $\in \CN^\d$ could be set according to the problem at hand or learned as hyperparameters. Also, note that kernel in \EQ{expkernel} is a particularization of \EQ{kerfil}, by setting $v\rrjj=0$.


\begin{figure}[tb!]
\begin{center}
\includegraphics[width=6cm, draft=false]{Fig1}
%\setlength\figureheight{5.5cm}
%\setlength\figurewidth{7cm}
%\input{figures/mifigura03.tex}
\end{center}
\vspace*{-.6cm}
\caption{Convolution system model to design a kernel for proper Gaussian processes.}
%\LABFIG{fig1}
%load prueba_fer2_30_50e6_3_1_50
%\vspace*{-.3cm}
\end{figure}







%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Hyperparameters computation}\LABSEC{ML}

%\subsection{The complex marginal likelihood}

In order to turn complex GPR into a powerful practical tool we must address the model selection problem, in the sense of setting the hyperparameters. A major advantage of GPR is that the hyperparameters can be estimated by maximizing the marginal likelihood \cite{Rasmussen06}. %A major advantage of the marginal likelihood in model selection is that it automatically incorporates a trade-off between model fitting and model complexity \cite{Rasmussen06}. %This is the reason why the marginal likelihood is valuable in solving the model selection problem. 
The marginal likelihood in (\ref{eq:marlikelihood}) yields a zero-mean complex proper Gaussian random vector with covariance matrix $\C$, and the log marginal likelihood is
\begin{equationw} \label{eq:CMML}
L=\log p\left(\yv_\n| \X_\n \right)=-\yv_\n\her\C\inv\yv_\n-\log\det\C-{n}\log \pi.
\end{equationw} 
In our model, $\C=\C(\boldsymbol{\theta})$ can be parameterized in terms of a set of hyperparameters, ${\theta}_{i}$, for the kernel $k(\x, \x')$ and the variance of the additive noise, $\sigma_\epsilon^2$, that can be treated as another hyperparameter. 
%Therefore, matrix $\C(\boldsymbol{\theta})$ and hence the log marginal likelihood $L\left(\boldsymbol{\theta}\right)$ can be explicitly written conditioned on the hyperparameters $\boldsymbol{\theta}$. 
 %Maximization of $L\left(\boldsymbol{\theta}\right)$ allows the optimization of those hyperparameters. 
Note that matrix $\C(\boldsymbol{\theta})$ is a covariance matrix, so it is Hermitian and $L\left(\boldsymbol{\theta}\right)$ is a real function. In the maximization of \EQ{CMML} its derivative is needed. However, the Cauchy-Riemann conditions dictate that nonconstant real-valued functions that are defined on complex domains are not holomorphic, so the complex derivative cannot be used. Instead, we must seek formal or {\it Wirtinger derivatives} \cite{Bouboulis11}, \cite{Hjorungnes07}. Also, since $\C(\boldsymbol{\theta})$ is Hermitian, it has a structure or pattern, therefore generalized complex-valued matrix derivatives need to be considered as follows \cite{Hjorungnes08}. Define %the following function
 \begin{equation} \label{eq:CMML2}
g(\hat\C,\hat\C^{*})=-\yv_\n\her\hat\C\inv\yv_\n-\log\det\hat\C-{n}\log \pi,
\end{equation} 
where $\hat\C$ is a matrix with independent components, i.e., not Hermitian. It is immediate to see that if $\hat\C$ is substituted with $\C$ in $g(\hat\C,\hat\C^{*})$, the marginal likelihood is produced, i.e.,
\begin{equation}
L\left(\boldsymbol{\theta}\right)=\left.g(\hat\C,\hat\C^{*})\right|_{\hat\C=\C(\boldsymbol{\theta})}.
\end{equation} 
Now, since $\hat\C$ is unpatterned, applying the chain rule leads to 
\begin{align} 
\frac{\partial {L}}{\partial \theta_i} &=\left.\frac{\partial {g(\hat\C,\hat\C^{*})}}{\partial \theta_i}\right|_{\hat\C=\C(\boldsymbol{\theta})} \nonumber\\
&=\sum_{k=1}^{n}{\sum_{l=1}^{n}\left(\left.{\frac{\partial {g(\hat\C,\hat\C^{*})}}{\partial (\hat{\C})_{lk}}}\right|_{\hat\C=\C(\boldsymbol{\theta})}\hspace{-0.2cm}{\frac{\partial {({\C})_{lk}}}{\partial \theta_i}} \right)}\nonumber\\
&={\mathrm{Tr}\left(\left(\left.{\frac{\partial {g(\hat\C,\hat\C^{*})}}{\partial \hat\C}}\right|_{\hat\C=\C(\boldsymbol{\theta})}\right)^\top{\frac{\partial {\C}}{\partial \theta_i}}\right)}.\label{eq:derivad}
\end{align}
%where $(\hat{\C})_{lk}$ and $({\C})_{lk}$ are the $(l,k)$ entries of matrices $\hat\C$ and $\C$, respectively.

The derivative of the first term of $g(\hat\C,\hat\C^{*})$ in (\ref{eq:CMML2}) yields
\begin{align} \label{eq:1term}
\frac{\partial {}}{\partial \hat\C}\left(-\yv_\n\her\hat\C\inv\yv_\n\right)&=-\frac{\partial {}}{\partial \hat\C}\mathrm{Tr}\left(-\yv_\n\her\hat\C\inv\yv_\n\right)\nonumber\\
%&=&-\frac{\partial {}}{\partial\C}\mathrm{Tr}\left(\yv\yv\her\C\inv\right)\nonumber\\
&=(\hat\C^{\top})^{-1}\left(\yv_\n\yv_\n\her\right)^\top(\hat\C^{\top})^{-1}.
\end{align}
The derivative of the second term of of $g(\hat\C,\hat\C^{*})$ in (\ref{eq:CMML2}) yields
\begin{equation} \label{eq:2term}
\frac{\partial {}}{\partial \hat\C}\left(-\log\det\hat\C\right)=-(\hat\C^{\top})^{-1}.
\end{equation}
Substitution of (\ref{eq:1term}) and (\ref{eq:2term}) in (\ref{eq:derivad}) yield
\begin{align} \label{eq:gradientefinal}
\frac{\partial {L}}{\partial \theta_i} ={\mathrm{Tr}\left(\left(\C^{-1}\yv_\n\yv_\n\her\C^{-1}-\C^{-1}\right){\frac{\partial {\C}}{\partial \theta_i}}\right)},
%\nonumber\\
%&={\mathrm{Tr}\left[\left(\C\inv\left(\yv_\n\yv_\n\her\right)\C\inv-\C\inv\right){\frac{\partial {\C}}{\partial \theta_i}}\right]}.
\end{align}
where now the pattern of matrix $\C$ does not yield further simplifications. 
The term ${\partial {\C}}/{\partial \theta_i}$ depends on the chosen covariance function.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{experiments}\LABSEC{Exp}

\subsection{Non-null Cross-covariances}

We propose the following example where we randomly sampled a proper complex GP, to later learn it. We generated samples for both the real and the imaginary parts in $[-6,5]$. The GP covariance matrix used was the one in \EQ{kerfil}, constructed from the corresponding filters in \EQ{fil}. The kernel hyperparameters were set as $\gamma=1.125$ and $\mu=2+2\j$. Values $v\rr$ and $v\rrjj$ were set to 1. The real part of the sample function is shown in \FIG{fig1} (top), while the imaginary part is included in \FIG{fig2} (top). Circular complex Gaussian noise with $\sigma_{\epsilon}=0.1$ was added to represent measurement uncertainty and $\n=200$ training noisy samples were randomly chosen. These samples has been depicted as circles in \FIG{fig1} and \FIG{fig2} (top). 
%We generated a set of $6400$ complex input points $\{x(i)| i = 1, . . . , 6400\}$ by taking $80$ samples for both the real and the imaginary parts linearly spaced in $[-6,4]$. Then we calculated the corresponding filters as in \EQ{fil} to yield the kernel in \EQ{kerfil}. The kernel hyperparameters were set as $\gamma=1.125$, $v\rrjj=v\rrjj=1$ and $\mu=2+2\j$. By stimulating the filters with zero-mean unit variance Gaussian noises we constructed a proper complex GP. Note that both the real and imaginary parts of the drawn process are 2-dimensional surfaces over the real and the imaginary parts of the inputs. The real part is shown in \FIG{fig1} (top), while the imaginary part is shown in \FIG{fig2} (top). Circular complex noise with $\sigma_{\epsilon}=0.1$ was added to represent measurement uncertainty and $\n=200$ training noisy samples were randomly chosen. 
Maximization of the log marginal likelihood in (\ref{eq:CMML}) using (\ref{eq:gradientefinal}) yielded the following estimated values of the hyperparameters: $\hat\gamma=1.0857$, $\hat\mu=1.9671+\j1.8690$, and $\hat\sigma_{\epsilon}=0.0984$. Then, we found the mean (\ref{eq:mediaCPGP}) and variance (\ref{eq:varCPGP}) of the predictive distribution using those training samples and the estimated values of the hyperparameters. The real and the imaginary parts of the predictive mean (\ref{eq:mediaCPGP}) are depicted in \FIG{fig1} (bottom) and \FIG{fig2} (bottom), respectively. The mean squared error of the estimation was $10\log_{10}MSE =-12.52$ dB, computed for 6400 inputs. 

To further study the accuracy of the prediction, we include in \FIG{fig3} and \FIG{fig4} two slices of the surfaces in \FIG{fig1} and \FIG{fig2}, respectively. In \FIG{fig3} the real part of the sample function of the process is plotted (dashed line) versus the real part of the input, the imaginary part of the input was fixed to the value $\Im(x)=4.4430$. The real part of the prediction in (\ref{eq:mediaCPGP}) is depicted in thick red line, along with the grey shaded area that represents the pointwise mean plus and minus two times the standard deviation. The blue circles mark the training samples. Similarly, in \FIG{fig4} the imaginary part of the input was fixed to the value $\Im(x)=-0.5696$ and the imaginary part of the sample function of the process is plotted (dashed line) versus the real part of the input. Again, the imaginary part of the prediction and the training samples are depicted. Four instances of the posterior (\ref{eq:posteriorf}) are plotted in thin blue lines in Figs. \ref{fig:fig3} and \ref{fig:fig4}.
\begin{figure}[tb!]
\begin{center}
\includegraphics[width=8.6cm, draft=false]{2Dpred_real}
%\setlength\figureheight{5.5cm}
%\setlength\figurewidth{7cm}
%\input{figures/mifigura03.tex}
\end{center}
\vspace*{-.6cm}
\caption{Real part of the output $y(\x)$ (top) and real part of the mean estimation (\ref{eq:mediaCPGP}) (bottom) versus the real and imaginary parts of the input, $\Re(x)$ and $\Im(x)$, respectively. The training samples are depicted as blue circles.}
\LABFIG{fig1}
%load prueba_fer2_30_50e6_3_1_50
%\vspace*{-.3cm}
\end{figure}

%
\begin{figure}[tb!]
\begin{center}
\includegraphics[width=8.6cm, draft=false]{2Dpred_imag}
%\setlength\figureheight{5.5cm}
%\setlength\figurewidth{7cm}
%\input{figures/mifigura03.tex}
\end{center}
\vspace*{-.6cm}
\caption{Imaginary part of the output $y(\x)$ (top) and imaginary part of the mean estimation (\ref{eq:mediaCPGP}) (bottom) versus the real and imaginary parts of the input, $\Re(x)$ and $\Im(x)$, respectively. The training samples are depicted as blue circles.}
\LABFIG{fig2}
%load prueba_fer2_30_50e6_3_1_50
%\vspace*{-.3cm}
\end{figure}
%
\begin{figure}[tb!]
\begin{center}
\includegraphics[width=8.5cm, draft=false]{Fig4}
%\setlength\figureheight{5cm}
%\setlength\figurewidth{7cm}
%\input{figures/1Dpred_re.tex}
\end{center}
\vspace*{-.6cm}
\caption{Real part of the output (thick black dashed line) and real part of the predictive mean (\ref{eq:mediaCPGP}) (thick red line) versus $\Re(x)$ for $\Im(x)=4.4430$. Training samples are depicted as blue circles. Four instances of the posterior (\ref{eq:posteriorf}) are plotted in thin blue lines.}
%load prueba_fer2_30_50e6_3_1_50
%\vspace*{-.3cm}
\LABFIG{fig3}
\end{figure}
%
\begin{figure}[tb!]
\begin{center}
\includegraphics[width=8.5cm, draft=false]{Fig5}
%\setlength\figureheight{5cm}
%\setlength\figurewidth{7cm}
%\input{figures/1Dpred_im.tex}
\end{center}
\vspace*{-.6cm}
\caption{Imaginary part of the output (thick black dashed line) and imaginary part of the predictive mean (\ref{eq:mediaCPGP}) (thick red line) versus $\Re(x)$ for $\Im(x)=-0.5696$. Training samples are depicted as blue circles. Four instances of the posterior (\ref{eq:posteriorf}) are plotted in thin blue lines.}
%load prueba_fer2_30_50e6_3_1_50
%\vspace*{-.3cm}
\LABFIG{fig4}
\end{figure}

%$k_1$$U_1$ that once convolved provide the real and imaginary parts of the outputs. 
%\begin{align}
%       
%\end{align}
%
%\k_{G}(\x,\x')

\subsection{Null Cross-covariances}
The performance of the proposed proper complex GPR was tested in the context of a nonlinear channel equalization task. Two nonlinear channels were considered, as in \cite{Bouboulis11} and \cite{Bouboulis12}. Both channels consisted of a linear filter $t(n)=(-0.9+0.8\j)\cdot s(n)+(0.6-0.7\j)\cdot s(n-1)$ and a memoryless nonlinearity. The nonlinearity was $q(n)=t(n)+(0.1+0.15\j)\cdot t^{2}(n)+ (0.06+0.05\j)\cdot t^{3}(n)$ for the first case (labeled as {\it soft nonlinear channel}), and $q(n)=t(n)+(0.2+0.25\j)\cdot t^{2}(n)+ (0.12+0.09\j)\cdot t^{3}(n)$ for the second case (labeled as {\it strong nonlinear channel}). The input signals had the form $s(n)=0.70(\sqrt{1-\rho^{2}}X(n)+\j\rho Y(n))$, as in \cite{Bouboulis11} and \cite{Bouboulis12}, and $X(n)$ and $Y(n)$ were Gaussian random variables. Note that the real and the imaginary parts of the input signals were generated independently and, therefore, had null cross-covariances. Also note that the input signals are circular for $\rho=1/\sqrt{2}$ and highly noncircular if $\rho$ approaches $0$ or $1$. At the receiver end of the channel, the signal $q(n)$ was corrupted by additive white circular Gaussian noise with the SNR set to 16 dB.

The aim of a channel equalization task is to construct an inverse filter, which acts on the received signal $r(t)$ and reproduces the original input signal $s(n)$ as close as possible. To this end, the inputs to the equalizer were the sets of samples $\vect{x}=[r(n+D),r(n+D-1),\cdots,r(n+D-L+1)]^{\top}$, where $L>0$ is the filter length and $D$ is the equalization time delay. %These parameters were set to $L=5$ and $D=2$, as in 

Experiments were conducted as in \cite{Bouboulis11} and \cite{Bouboulis12}, where $L=5$ and $D=2$, on a set of 3000 samples of the input signal considering both the circular and the noncircular ($\rho=0.1$) cases and the described nonlinear channels ({soft} and {strong}). In all cases the results were averaged over 100 trials where the input signals $s(n)$ and noise output were generated randomly. In Figs. \ref{fig:fig5} to \ref{fig:fig8} we include the MSE along training samples for several methods and the channels and inputs described. As in \cite{Bouboulis11} and \cite{Bouboulis12} the predicted outputs for the training inputs were used to compute the error. The MSE value depicted for each sample is the averaged MSE for the last $100$ samples as in \cite{boubouliscode}\footnote{In the code the average was performed for the last 500 samples.}. %And the MSE shown is 
%show the good performance in all the cases, soft or strong nonlinear channels and circular or noncircular signals. 
%
For the sake of comparison we depict the results for the NCLMS, the NCKLMS2 in \cite{Bouboulis11} and the ACKLMS \cite{Bouboulis12} algorithms, that use the complex Gaussian kernel in \EQ{RK}. We used the code available in \cite{boubouliscode} to run these algorithms. For the ACKLMS algorithm all the parameters were set to the values described in \cite{Bouboulis12} ($\gamma=10^{2}$ and the step update parameter $1/8$). Also, for the NCKLMS2 algorithm and the soft nonlinear channel case the parameters were set to the values described in \cite{Bouboulis12} ($\gamma=10^{2}$ and the step update parameter $1/8$), whereas for the strong nonlinear channel case the parameters were set to the values described in \cite{Bouboulis11} ($\gamma=5^{2}$ and the step update parameter $1/4$). The novelty criterion was used for the sparsification of both the NCKLMS2 and the ACKLMS algorithms with $\delta_{1}=0.15$ and $\delta_{2}=0.2$. %Figs. \ref{fig:fig5} to \ref{fig:fig8} show the learning curves. 
 In Figs. \ref{fig:fig7} and \ref{fig:fig8} we observe stability problems in the learning process of both NCKLMS2 and ACKLMS, due to the kernel used. This problem was alleviated by using sparsification. We also used the independent kernel (\ref{eq:indepCGK}) with $\kappa_{\Rext}$ being the real Gaussian kernel in the NCKLMS2 approach, as proposed in \cite{Tobar12}. We labeled this algorithm as NCKLMS2-i. The tunable parameter was set to $\gamma=5^{2}$ and the step update parameter to $1/8$. The novelty criterion was again used for the sparsification. 

In this paper we propose the mean (\ref{eq:mediaCPGP}) of the predictive distribution for the proper complex GPR as the equalizer output. We used the kernel in \EQ{complexCV} where $k\rrrr(\x,\x')$$=k\jjjj(\x,\x')$$=\k_{G}(\x,\x')/2$ with $\k_{G}(\x,\x')$ as in (\ref{eq:expkernel}), and $k\rrjj(\x,\x')=0$, since the real and the imaginary parts of $s(n)$ had null cross-covariances. Two approaches were simulated. In one procedure, which we labeled as opt-CGPR, $1000$ input-output samples were randomly chosen to tune the kernel hyperparameter $\gamma$ so that the best possible results were obtained. The noise variance $\sigma_{\epsilon}^{2}$ in the kernel was set to the valued used in the experiments. Then, with those hyperparameter values, for each new input $\vect{x}$ the equalizer predicted the corresponding $s(n)$ as in (\ref{eq:mediaCPGP}) using all the previous input-output pairs as training set. It can be observed through the figures the remarkable good results of this proposal in all the cases, soft or strong nonlinear channels and circular or noncircular signals. Also, note that the MSE increases with the number of samples until reaches steady estate, as CGPR generalizes for an increasing number of input-output pairs.
In a second approach, labeled as CGPR, we tried to check the hyperparameter estimation capabilities of the proper complex GPR by means of the maximization of the log marginal likelihood (\ref{eq:CMML}) using (\ref{eq:gradientefinal}). The first $250$ samples were used as training set to estimate the hyperparameters ($\gamma$ and $\sigma_{\epsilon}$). Then, with those values for $\gamma$ and $\sigma_{\epsilon}^{2}$ fixed, and for each new input $\vect{x}$ the equalizer predicted the corresponding $s(n)$ as in (\ref{eq:mediaCPGP}) using all the previous input-output pairs as training set. The learning curves for this procedure are depicted in Figs. \ref{fig:fig5} to \ref{fig:fig8}. The results for this procedure, CGPR, not far from that of the opt-CGPR, considerably outperform the NCKLMS2 and ACKLMS and illustrate the good estimation of the hyperparameters. %Figs. \ref{fig:fig5} to \ref{fig:fig8} show that the proposed equalizers, CGPR and opt-CGPR, considerably outperform the NCKLMS2 and ACKLMS, and have the additional advantage of hyperparameters optimization capability. 

Finally, we proposed to improve the NCKLMS2 by using the kernels derived in \SEC{Cov}, in particular the kernel in (\ref{eq:expkernel}) with null imaginary part. This algorithm is labeled as NCKLMS2-G. We set all the parameters for this algorithm to the same values indicated in \cite{Bouboulis11} ($\gamma=5^{2}$ and the step update parameter $1/4$), and the novelty criterion was again used for the sparsification. While the NCKLMS2-i algorithm exhibits improved performance compared to the NCKLMS2 or ACKLMS algorithms, the NCKLMS2-G algorithm outperforms all of them with a smooth learning curve under all conditions. This proves the proposed complex kernel design procedure to be a valuable one. In this equalization problem the cross-covariance between real and imaginary parts of the signals to-be-learned, $s(n)$, is null. Therefore, we set the imaginary part of the kernel to zero with much better results. Also, the measure of similarity of inputs through the simple norm of the complex difference between inputs endows the kernel with useful properties such as isotropy and stationarity, better fitting the underlying model. % and resort to the ability of the complex GPR to learn the outputs the proposed kernel, which does not have imaginary part, is more suitable than those that have imaginary part as the complex Gaussian kernel \EQ{RK} or the independent kernel (\ref{eq:indepCGK}).
\begin{figure}[tb!]
\begin{center}
\includegraphics[width=8.6cm, draft=false]{soft_circ2.pdf}
\end{center}
%\vspace*{-.6cm}
\caption{Learning curves for NKCLMS2, NCKLMS2-i, NKCLMS2-G, AKCLMS, the proper complex GPR (CGPR) and the optimized CGPR (opt-CGPR) for the soft nonlinear channel equalization problem for the circular input case.}
\LABFIG{fig5}
%load prueba_fer2_30_50e6_3_1_50
%\vspace*{-.3cm}
\end{figure}
\begin{figure}[tb!]
\begin{center}
\includegraphics[width=8.6cm, draft=false]{circular2.pdf}
\end{center}
%\vspace*{-.6cm}
\caption{Learning curves for NKCLMS2, NCKLMS2-i, NKCLMS2-G, AKCLMS, the proper complex GPR (CGPR) and the optimized CGPR (opt-CGPR) for the strong nonlinear channel equalization problem for the circular input case.}
\LABFIG{fig6}
%load prueba_fer2_30_50e6_3_1_50
%\vspace*{-.3cm}
\end{figure}

\begin{figure}[tb!]
\begin{center}
\includegraphics[width=8.6cm, draft=false]{soft_nocirc2.pdf}
%\setlength\figureheight{5.5cm}
%\setlength\figurewidth{7cm}
%\input{figures/mifigura03.tex}
\end{center}
\vspace*{-.6cm}
\caption{Learning curves for NKCLMS2, NCKLMS2-i, NKCLMS2-G, AKCLMS, the proper complex GPR (CGPR) and the optimized CGPR (opt-CGPR) for the soft nonlinear channel equalization problem for the noncircular input case ($\rho=0.1$).}
\LABFIG{fig7}
%load prueba_fer2_30_50e6_3_1_50
%\vspace*{-.3cm}
\end{figure}
\begin{figure}[tb!]
\begin{center}
\includegraphics[width=8.6cm, draft=false]{nocircular2.pdf}
%\setlength\figureheight{5.5cm}
%\setlength\figurewidth{7cm}
%\input{figures/mifigura03.tex}
\end{center}
\vspace*{-.6cm}
\caption{Learning curves for NKCLMS2, NCKLMS2-i, NKCLMS2-G, AKCLMS, the proper complex GPR (CGPR) and the optimized CGPR (opt-CGPR) for the strong nonlinear channel equalization problem for the noncircular input case ($\rho=0.1$).}
\LABFIG{fig8}
%load prueba_fer2_30_50e6_3_1_50
%\vspace*{-.3cm}
\end{figure}


\section{Conclusions}

In this paper we prove that the kernels for complex-valued regression proposed in the literature correspond to \emph{proper} complex outputs. Since proper complex-valued Gaussian processes have been widely studied, we propose to resort to the Gaussian processes for regression to study this problem. On the one hand, we end with a new solution that it is endowed with the advantages of GPs: a probabilistic output and the availability of estimating the hyperparameters by optimization. On the other, we derive an expression for the reproducing kernel in the proper complex case. We conclude that it is complex valued if the cross-covariance of the real and imaginary parts of the outputs is not null. The imaginary part of the covariance matrix is skew-symmetric, and the full matrix must be semidefinite positive, but the kernel used in the real and imaginary parts need not to be the same. Besides, if the cross-covariance cancels, the kernel is real-valued. The MOL point of view bears this out. Therefore, we prove the kernel to exhibit a much more flexible structure than proposed in previous works. We end designing a kernel where the use of an appropriate metric yields an isotropic and stationary kernel. The proposed methods developed by fully using these key results exhibit a remarkable good performance compared to previous ones. 

%show that a complex-valued formulation of the GPR can be derived for the proper case. As a result of this statistical description, we discuss the structure that a complex covariance matrix should have for the proper case and how kernels can be derived from real kernels. Also, we develop the useful feature of the GPR to learn hyperparameters from the optimization of the marginal likelihood to be also available in the proper complex GPR. These results make proper CGPR a more powerful approach compared to previous proposals for complex values, based on kernel least square, where the hyperparameters selection and the kernel design were somehow undefined. The examples included illustrate the good behavior of the proper CGPR and the performance in the hyperparameter estimation when learning a proper Gaussian process and in the context of nonlinear channel equalization.



% if have a single appendix:
\appendix[Covariance Function]
We follow here a procedure similar to that in \cite{Boyle05}. Consider two independent, real, stationary, Gaussian white noise processes $S\rr(\x)$ and $S\jj(\x)$, where $\x\in\CN^\d$, producing an output $Y(\x)=U(\x)+W(\x)$, where $W(\x)$ is a stationary Gaussian white noise, and $U(\x)$ is defined by the sum of convolutions
\begin{align}
U(\x)&=\left(h\rr(\x)+\j h\rrjj(\x)\right) \star S\rr(\x)+\left(h\rr(\x)+\j h\jjrr(\x)\right) \star S\jj(\x)\nonumber\\
&=\sum_{m=1}^4 h_m(\x) \star S_m(\x),
\end{align}
where we have used the following notation: $h_1(\x)=h\rr(\x)$, $h_2(\x)=\j h\rrjj(\x)$, $h_3(\x)=h\rr(\x)$, $h_4(\x)=\j h\jjrr(\x)$, $S_1(\x)=S_2(\x)= S\rr(\x)$, and $S_3(\x)=S_4(\x)= S\jj(\x)$.
The covariance of $Y(\x)$ is derived as follows:
\begin{align}
\C(\x_a,\x_b)=\C_{U}(\x_a,\x_b)+\sigma_W^2\delta_{ab},
\end{align}
where $\sigma_W^2$ is the variance of $W(\x)$, and
\begin{align} \LABEQ{app2}
\C_{U}(\x_a,\x_b)&=\E\left[U\rr(\x_a)U\rr^{*}(\x_b)\right]\nonumber\\
&=\E\left[\sum_{m=1}^4\int_{\CN^\d}{h_m(\mathbf{\alpha})S_m(\x_a-\mathbf{\alpha})\dfd{\d}{\mathbf{\alpha}}}\right.\nonumber\\&\cdot\left.\sum_{n=1}^4\int_{\CN^\d}{h_n^{*}(\mathbf{\beta})S_n(\x_b-\mathbf{\beta})\dfd{\d}{\mathbf{\beta}}}\right]
\nonumber\\
&=\sum_{m=1}^4\sum_{n=1}^4\left\{\int_{\CN^\d}\int_{\CN^\d}h_m(\mathbf{\alpha})h_n^{*}(\mathbf{\beta})\right.\nonumber\\&\cdot\left.\E\left[S_m(\x_a-\mathbf{\alpha})S_n(\x_b-\mathbf{\beta})\right]\dfd{\d}{\mathbf{\alpha}}\dfd{\d}{\mathbf{\beta}}\right\}.
\end{align}
Since $S_1(\x)=S_2(\x)= S\rr(\x)$, and $S_3(\x)=S_4(\x)= S\jj(\x)$, processes $S_m(\x_a-\mathbf{\alpha})$ and $S_n(\x_b-\mathbf{\beta})$ covary only if $m,n\in\{1,2\}$ or $m,n\in\{3,4\}$, and $(\x_a-\mathbf{\alpha})=(\x_b-\mathbf{\beta})$. In such cases, $\E\left[S_m(\x_a-\mathbf{\alpha})S_n(\x_b-\mathbf{\beta})\right]=\delta(\mathbf{\alpha}-(\x_a-\x_b+\mathbf{\beta}))$, and the integrals in \EQ{app2} yield
\begin{align}
&\int_{\CN^\d}\int_{\CN^\d}h_m(\mathbf{\alpha})h_n^{*}(\mathbf{\beta})\delta(\mathbf{\alpha}-(\x_a-\x_b+\mathbf{\beta}))\dfd{\d}{\mathbf{\alpha}}\dfd{\d}{\mathbf{\beta}}\nonumber\\
&=\int_{\CN^\d}h_m(\mathbf{\beta}+(\x_a-\x_b))h_n^{*}(\mathbf{\beta})\dfd{\d}{\mathbf{\beta}}\nonumber\\
&=\int_{\CN^\d}h_m(\mathbf{\beta}+\mathbf{d}_\x)h_n^{*}(\mathbf{\beta})\dfd{\d}{\mathbf{\beta}},
\end{align}
where $\mathbf{d}_\x=\x_a-\x_b$.
Hence,
\begin{align} \LABEQ{app1}
\C_{U}(\x_a,\x_b)&=\sum_{m=1}^2\sum_{n=1}^2\int_{\CN^\d}h_m(\mathbf{\beta}+\mathbf{d}_\x)h_n^{*}(\mathbf{\beta})\dfd{\d}{\mathbf{\beta}}\nonumber\\
&+\sum_{m=3}^4\sum_{n=3}^4\int_{\CN^\d}h_m(\mathbf{\beta}+\mathbf{d}_\x)h_n^{*}(\mathbf{\beta})\dfd{\d}{\mathbf{\beta}}.
\end{align}
If we set the kernels to parameterized exponentials as $h_i(\x)=(\j)^{c_{i}}v_{i}\exp(-(\x-\mu_i)\her(\x-\mu_i)/\gamma)$, integrals in \EQ{app1} are as follows
 \begin{align}\LABEQ{app3}
&\int_{\CN^\d}h_m(\mathbf{\beta}+\mathbf{d}_\x)h_n^{*}(\mathbf{\beta})\df{\mathbf{\beta}^\d}\nonumber
\\=&(\j)^{c_{m}}(-\j)^{c_{n}}v_{m}v_{n}\int_{\CN^\d}\exp\left(-\frac{(\mathbf{\beta}-\mu_n)\her(\mathbf{\beta}-\mu_n)}{\gamma}\right)\nonumber\\\cdot&\exp\left(-\frac{(\mathbf{\beta}+\mathbf{d}_\x-\mu_m)\her(\mathbf{\beta}+\mathbf{d}_\x-\mu_m)}{\gamma}\right)\df{\mathbf{\beta}^\d}\nonumber\\
=&(\j)^{c_{m}}(-\j)^{c_{n}}v_{m}v_{n}\left(\int_{\CN^\d}\exp\left(-\frac{(\mathbf{\beta}-\hat\beta)\her(\mathbf{\beta}-\hat\beta)}{0.5\gamma}\right)\df{\mathbf{\beta}^\d}\right)\nonumber\\&\cdot\exp\left(-\frac{(\mathbf{d}_\x-\mu_m+\mu_n)\her(\mathbf{d}_\x-\mu_m+\mu_n)}{2\gamma}\right)\nonumber\\
=&(\j)^{c_{m}}(-\j)^{c_{n}}\left(\frac{\pi\gamma}{2}\right)^{\d}v_{m}v_{n}\nonumber\\\cdot&\exp\left(-\frac{(\mathbf{d}_\x-\mu_m+\mu_n)\her(\mathbf{d}_\x-\mu_m+\mu_n)}{2\gamma}\right),
\end{align}
where $\hat\beta=0.5(\mu_n-(\mathbf{d}_\x-\mu_m))$. 

We propose the following parameter values. For $h_1(\x)=h_3(\x)=h\rr(\x)$, we set $v_{1}=v_{3}=v\rr$, $\mu_{1}=\mu_{3}=0$ and $c_{1}=c_{3}=0$. For $h_2(\x)=\j h\rrjj(\x)$, we set $\mu_{2}=\mu$, $v_{2}=v\rrjj$ and $c_{2}=1$. And for $h_4(\x)=\j h\jjrr(\x)$, we set $\mu_{4}=-\mu$, $v_{4}=-v\rrjj$ and $c_{4}=1$. By making use of this values and \EQ{app3}, \EQ{app1} yields
\begin{align} 
&\C_{U}(\x_a,\x_b)=\C_{U}(\mathbf{d}_\x)\nonumber\\
%&\int_{\CN^\d}h_{1}(\mathbf{\beta}+(\mathbf{d}_\x))h_{1}(\mathbf{\beta})\df{\mathbf{\beta}^\d}+\int_{\CN^\d}h_{1}(\mathbf{\beta}+(\mathbf{d}_\x))h_{2}(\mathbf{\beta})\df{\mathbf{\beta}^\d}\nonumber\\
%&+\int_{\CN^\d}h_{2}(\mathbf{\beta}+\mathbf{d}_\x)h_{1}(\mathbf{\beta})\df{\mathbf{\beta}^\d}+\int_{\CN^\d}h_{2}(\mathbf{\beta}+\mathbf{d}_\x)h_{2}(\mathbf{\beta})\df{\mathbf{\beta}^\d}\nonumber\\
%&+\int_{\CN^\d}h_{3}(\mathbf{\beta}+\mathbf{d}_\x)h_{3}(\mathbf{\beta})\df{\mathbf{\beta}^\d}+\int_{\CN^\d}h_{3}(\mathbf{\beta}+\mathbf{d}_\x)h_{4}(\mathbf{\beta})\df{\mathbf{\beta}^\d}\nonumber\\
%&+\int_{\CN^\d}h_{4}(\mathbf{\beta}+\mathbf{d}_\x)h_{3}(\mathbf{\beta})\df{\mathbf{\beta}^\d}+\int_{\CN^\d}h_{4}(\mathbf{\beta}+\mathbf{d}_\x)h_{4}(\mathbf{\beta})\df{\mathbf{\beta}^\d}\nonumber\\
&=\left(\frac{\pi\gamma}{2}\right)^{\d}v\rr^{2}\exp\left(-\frac{\mathbf{d}_\x\her\mathbf{d}_\x}{2\gamma}\right)\nonumber\\
&-\j\left(\frac{\pi\gamma}{2}\right)^{\d}v\rr v\rrjj\exp\left(-\frac{(\mathbf{d}_\x+\mu)\her(\mathbf{d}_\x+\mu)}{2\gamma}\right)\nonumber\\
&+\j\left(\frac{\pi\gamma}{2}\right)^{\d}v\rr v\rrjj\exp\left(-\frac{(\mathbf{d}_\x-\mu)\her(\mathbf{d}_\x-\mu)}{2\gamma}\right)\nonumber\\
&+\left(\frac{\pi\gamma}{2}\right)^{\d}v\rrjj^{2}\exp\left(-\frac{\mathbf{d}_\x\her\mathbf{d}_\x}{2\gamma}\right)+\left(\frac{\pi\gamma}{2}\right)^{\d}v\rr^{2}\exp\left(-\frac{\mathbf{d}_\x\her\mathbf{d}_\x}{2\gamma}\right)\nonumber\\
&+\j\left(\frac{\pi\gamma}{2}\right)^{\d}v\rr v\rrjj\exp\left(-\frac{(\mathbf{d}_\x-\mu)\her(\mathbf{d}_\x-\mu)}{2\gamma}\right)\nonumber\\
&-\j\left(\frac{\pi\gamma}{2}\right)^{\d}v\rr v\rrjj\exp\left(-\frac{(\mathbf{d}_\x+\mu)\her(\mathbf{d}_\x+\mu)}{2\gamma}\right)\nonumber\\
&+\left(\frac{\pi\gamma}{2}\right)^{\d}v\rrjj^{2}\exp\left(-\frac{\mathbf{d}_\x\her\mathbf{d}_\x}{2\gamma}\right)\nonumber\\
&=\left(\frac{\pi\gamma}{2}\right)^{\d}(2v\rr^{2}+2v\rrjj^{2})\exp\left(-\frac{\mathbf{d}_\x\her\mathbf{d}_\x}{2\gamma}\right)\nonumber\\
&+\j\left(\frac{\pi\gamma}{2}\right)^{\d}2v\rr v\rrjj\left(\exp\left(-\frac{(\mathbf{d}_\x-\mu)\her(\mathbf{d}_\x-\mu)}{2\gamma}\right)\right.\nonumber\\
&-\left.\exp\left(-\frac{(\mathbf{d}_\x+\mu)\her(\mathbf{d}_\x+\mu)}{2\gamma}\right)\right).
\end{align}




% if have a single appendix:
%\appendix[Proof of the Zonklar Equations]
% or
%\appendix  % for no appendix heading
% do not use \section anymore after \appendix, only \section*
% is possibly needed

% use appendices with more than one appendix
% then use \section to start each appendix
% you must declare a \section before using any
% \subsection or using \label (\appendices by itself
% starts a section numbered zero.)
%

% use section* for acknowledgement
%\section*{Acknowledgment}
%The authors would like to thank...


% Can use something like this to put references on a page
% by themselves when using endfloat and the captionsoff option.
\ifCLASSOPTIONcaptionsoff
  \newpage
\fi



% trigger a \newpage just before the given reference
% number - used to balance the columns on the last page
% adjust value as needed - may need to be readjusted if
% the document is modified later
%\IEEEtriggeratref{8}
% The "triggered" command can be changed if desired:
%\IEEEtriggercmd{\enlargethispage{-5in}}

% references section

% can use a bibliography generated by BibTeX as a .bbl file
% BibTeX documentation can be easily obtained at:
% http://www.ctan.org/tex-archive/biblio/bibtex/contrib/doc/
% The IEEEtran BibTeX style support page is at:
% http://www.michaelshell.org/tex/ieeetran/bibtex/
%\bibliographystyle{IEEEtran}
% argument is your BibTeX string definitions and bibliography database(s)
%\bibliography{IEEEabrv,../bib/paper}
%
% <OR> manually copy in the resultant .bbl file
% set second argument of \begin to the number of references
% (used to reserve space for the reference number labels box)

\bibliographystyle{IEEEtran}
\bibliography{IEEEabrv,CGPR,murilloGP,SSCDMA}

%\begin{thebibliography}{1}
%
%\bibitem{IEEEhowto:kopka}
%H.~Kopka and P.~W. Daly, \emph{A Guide to \LaTeX}, 3rd~ed.\hskip 1em plus
%  0.5em minus 0.4em\relax Harlow, England: Addison-Wesley, 1999.
%
%\end{thebibliography}

% biography section
% 
% If you have an EPS/PDF photo (graphicx package needed) extra braces are
% needed around the contents of the optional argument to biography to prevent
% the LaTeX parser from getting confused when it sees the complicated
% \includegraphics command within an optional argument. (You could create
% your own custom macro containing the \includegraphics command to make things
% simpler here.)
%\begin{IEEEbiography}[{\includegraphics[width=1in,height=1.25in,clip,keepaspectratio]{mshell}}]{Michael Shell}
% or if you just want to reserve a space for a photo:

%\begin{IEEEbiography}{Michael Shell}
%Biography text here.
%\end{IEEEbiography}

% if you will not have a photo at all:
%\begin{IEEEbiographynophoto}{John Doe}
%Biography text here.
%\end{IEEEbiographynophoto}

% insert where needed to balance the two columns on the last page with
% biographies
%\newpage

%\begin{IEEEbiographynophoto}{Jane Doe}
%Biography text here.
%\end{IEEEbiographynophoto}

% You can push biographies down or up by placing
% a \vfill before or after them. The appropriate
% use of \vfill depends on what kind of text is
% on the last page and whether or not the columns
% are being equalized.

%\vfill

% Can be used to pull up biographies so that the bottom of the last one
% is flush with the other column.
%\enlargethispage{-5in}



% that's all folks
\end{document}

arxiv2015.bbl
% Generated by IEEEtran.bst, version: 1.13 (2008/09/30)
\begin{thebibliography}{10}
\providecommand{\url}[1]{#1}
\csname url@samestyle\endcsname
\providecommand{\newblock}{\relax}
\providecommand{\bibinfo}[2]{#2}
\providecommand{\BIBentrySTDinterwordspacing}{\spaceskip=0pt\relax}
\providecommand{\BIBentryALTinterwordstretchfactor}{4}
\providecommand{\BIBentryALTinterwordspacing}{\spaceskip=\fontdimen2\font plus
\BIBentryALTinterwordstretchfactor\fontdimen3\font minus
  \fontdimen4\font\relax}
\providecommand{\BIBforeignlanguage}[2]{{%
\expandafter\ifx\csname l@#1\endcsname\relax
\typeout{** WARNING: IEEEtran.bst: No hyphenation pattern has been}%
\typeout{** loaded for the language `#1'. Using the pattern for}%
\typeout{** the default language instead.}%
\else
\language=\csname l@#1\endcsname
\fi
#2}}
\providecommand{\BIBdecl}{\relax}
\BIBdecl

\bibitem{Schreier06}
P.~J. Schreier and L.~L. Scharf, \emph{Statistical Signal Processing of
  Complex-Valued Data. The Theory of Improper and Noncircular Signals}.\hskip
  1em plus 0.5em minus 0.4em\relax Cambridge, UK: Cambridge University Press,
  2010.

\bibitem{Neeser93}
F.~Neeser and J.~Massey, ``Proper complex random processes with applications to
  information theory,'' \emph{IEEE Tr. on Information Theory}, vol.~39, no.~4,
  pp. 1293--1302, Jul 1993.

\bibitem{Tulay11}
T.~Adali, P.~J. Schreier, and L.~L. Scharf, ``Complex-valued signal processing:
  The proper way to deal with impropriety.'' \emph{IEEE Trans. Signal
  Processing}, pp. 5101--5125, 2011.

\bibitem{hirose13}
A.~Hirose, \emph{Complex-Valued Neural Networks: Advances and Applications},
  ser. {IEEE} Press Series on Computational Intelligence.\hskip 1em plus 0.5em
  minus 0.4em\relax Wiley, 2013.

\bibitem{Scholkopf02}
B.~Sch{\"o}lkopf and A.~Smola, \emph{Learning with Kernels: Support Vector
  Machines, Regularization, Optimization, and Beyond}, ser. Adaptive
  computation and machine learning.\hskip 1em plus 0.5em minus 0.4em\relax
  Cambridge, Massachusetts: MIT Press, 2002.

\bibitem{Steinwart06}
I.~Steinwart, D.~Hush, and C.~Scovel, ``An explicit description of the
  reproducing kernel {H}ilbert spaces of {G}aussian {RBF} kernels,'' \emph{IEEE
  Trans. Information Theory}, vol.~52, no.~10, pp. 4635--4643, Oct 2006.

\bibitem{OgunfunmiP11}
T.~Ogunfunmi and T.~K. Paul, ``On the complex kernel-based adaptive filter.''
  in \emph{IEEE Int. Symp. on Circuits and Systems, ISCAS}, 2011, pp.
  1263--1266.

\bibitem{Bouboulis12}
P.~Bouboulis, S.~Theodoridis, and M.~Mavroforakis, ``The augmented complex
  kernel {LMS},'' \emph{IEEE Trans. Signal Processing}, vol.~60, no.~9, pp.
  4962--4967, 2012.

\bibitem{Tobar12}
F.~Tobar, A.~Kuh, and D.~Mandic, ``A novel augmented complex valued kernel
  {LMS},'' in \emph{IEEE 7th Sensor Array and Multichannel Signal Processing
  Workshop, SAM}, 2012, pp. 473--476.

\bibitem{Papaioannou14}
A.~Papaioannou and S.~Zafeiriou, ``Principal component analysis with complex
  kernel: The widely linear model,'' \emph{IEEE Trans. Neural Networks and
  Learning Systems}, vol.~25, no.~9, pp. 1719--1726, Sept 2014.

\bibitem{Bouboulis11}
P.~Bouboulis and S.~Theodoridis, ``Extension of wirtinger's calculus to
  reproducing kernel {H}ilbert spaces and the complex kernel {LMS},''
  \emph{IEEE Trans. Signal Processing}, vol.~59, no.~3, pp. 964--978, 2011.

\bibitem{Boloix14}
R.~Boloix-Tortosa, F.~Payan-Somet, and J.~Murillo-Fuentes, ``Gaussian processes
  regressors for complex proper signals in digital communications,'' in
  \emph{IEEE 8th Sensor Array and Multichannel Signal Processing Workshop,
  SAM}, June 2014, pp. 137--140.

\bibitem{Rasmussen06}
C.~E. Rasmussen and C.~K.~I. Williams, \emph{Gaussian Processes for Machine
  Learning}.\hskip 1em plus 0.5em minus 0.4em\relax Cambridge, Massachusetts:
  MIT Press, 2006.

\bibitem{PerezCruz13gp}
F.~P\'erez-Cruz, S.~Van~Vaerenbergh, J.~Murillo-Fuentes, M.~Lazaro-Gredilla,
  and I.~Santamaria, ``Gaussian processes for nonlinear signal processing: An
  overview of recent advances,'' \emph{IEEE Signal Processing Magazine},
  vol.~30, no.~4, pp. 40--50, July 2013.

\bibitem{Miller69}
K.~Miller, ``Complex gaussian processes,'' \emph{SIAM Review}, vol.~11, no.~4,
  pp. 544--567, 1969.

\bibitem{Alvarez12}
M.~A. \'{A}lvarez, L.~Rosasco, and N.~D. Lawrence, ``Kernels for vector-valued
  functions: A review,'' \emph{Found. Trends Mach. Learn.}, vol.~4, no.~3, pp.
  195--266, Mar. 2012.

\bibitem{Taubock12}
G.~Taubo{\"o}ck, ``Complex-valued random vectors and channels: Entropy,
  divergence, and capacity,'' \emph{IEEE Trans. Information Theory}, vol.~58,
  no.~5, pp. 2729--2744, May 2012.

\bibitem{Boyle05}
P.~Boyle and M.~Frean, ``Dependent gaussian processes,'' in \emph{In Advances
  in Neural Information Processing Systems 17}.\hskip 1em plus 0.5em minus
  0.4em\relax MIT Press, 2005, pp. 217--224.

\bibitem{Calder07}
C.~Calder and N.~Cressie, ``Some topics in convolution-based spatial
  modeling,'' in \emph{Proceedings of the 56th Session of the International
  Statistics Institute}, Lisbon, Portugal, August 2007.

\bibitem{Ehrlich70}
L.~Ehrlich, ``Complex matrix inversion versus real,'' \emph{Communications of
  the ACM}, vol.~13, no.~9, pp. 561--562, 1970.

\bibitem{Golub96}
G.~H. Golub and C.~F. Van~Loan, \emph{Matrix Computations (3rd Ed.)}.\hskip 1em
  plus 0.5em minus 0.4em\relax Baltimore, MD, USA: Johns Hopkins University
  Press, 1996.

\bibitem{Hjorungnes08}
A.~Hjorungnes and D.~Palomar, ``Patterned complex-valued matrix derivatives,''
  in \emph{IEEE 5th Sensor Array and Multichannel Signal Processing Workshop,
  SAM 2008}, July 2008, pp. 293--297.

\bibitem{MurilloGP09}
J.~Murillo-Fuentes and F.~P\'erez-Cruz, ``Gaussian process regressors for
  multiuser detection in {DS-CDMA} systems,'' \emph{IEEE Trans.
  Communications}, vol.~57, no.~8, pp. 2339--2347, 2009.

\bibitem{PerezCruz08b}
F.~P\'erez-Cruz, J.~Murillo-Fuentes, and S.~Caro, ``Nonlinear channel
  equalization with {G}aussian {P}rocesses for {R}egression,'' \emph{IEEE
  Trans. Signal Processing}, vol.~56, no.~10, pp. 5283--5286, Oct. 2008.

\bibitem{Ollila08}
E.~Ollila, ``On the circularity of a complex random variable,'' \emph{IEEE
  Signal Processing Letters}, vol.~15, pp. 841--844, 2008.

\bibitem{PerezCruz08}
F.~Perez-Cruz, J.~Murillo-Fuentes, and S.~Caro, ``Nonlinear channel
  equalization with {G}aussian processes for regression,'' \emph{IEEE Trans.
  Signal Processing}, vol.~56, no.~10, pp. 5283--5286, Oct 2008.

\bibitem{Micchelli05}
C.~A. Micchelli and M.~Pontil, ``On learning vector-valued functions,''
  \emph{Neural Computation}, vol.~17, no.~1, pp. 177--204, 2005.

\bibitem{Hjorungnes07}
A.~Hjorungnes and D.~Gesbert, ``Complex-valued matrix differentiation:
  Techniques and key results,'' \emph{IEEE Trans. Signal Processing}, vol.~55,
  no.~6, pp. 2740--2746, June 2007.

\bibitem{boubouliscode}
P.~Bouboulis, ``{The Augmented Complex Kernel LMS Matlab and C Code},''
  \url{http://bouboulis.mysch.gr/kernels.html}, [Online].

\end{thebibliography}
