
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%       The Loss Rank Principle for Model Selection         %%
%%           Marcus Hutter, Start: September 2006            %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\documentclass[12pt,twoside]{article}
\usepackage{latexsym,amsmath,amssymb,hyperref}
\topmargin=-10mm  \oddsidemargin=5mm \evensidemargin=5mm
\textwidth=15cm \textheight=22cm
\sloppy\lineskip=0pt

%-------------------------------%
%       Macro-Definitions       %
%-------------------------------%
\def\eqvsp{}  \newdimen\paravsp  \paravsp=1.3ex

\def\,{\mskip 3mu} \def\>{\mskip 4mu plus 2mu minus 4mu} \def\;{\mskip 5mu plus 5mu} \def\!{\mskip-3mu}
\def\dispmuskip{\thinmuskip= 3mu plus 0mu minus 2mu \medmuskip=  4mu plus 2mu minus 2mu \thickmuskip=5mu plus 5mu minus 2mu}
\def\textmuskip{\thinmuskip= 0mu                    \medmuskip=  1mu plus 1mu minus 1mu \thickmuskip=2mu plus 3mu minus 1mu}
\textmuskip
\def\beq{\eqvsp\dispmuskip\begin{equation}}    \def\eeq{\eqvsp\end{equation}\textmuskip}
\def\beqn{\eqvsp\dispmuskip\begin{displaymath}}\def\eeqn{\eqvsp\end{displaymath}\textmuskip}
\def\bqa{\eqvsp\dispmuskip\begin{eqnarray}}    \def\eqa{\eqvsp\end{eqnarray}\textmuskip}
\def\bqan{\eqvsp\dispmuskip\begin{eqnarray*}}  \def\eqan{\eqvsp\end{eqnarray*}\textmuskip}

\newtheorem{theorem}{Theorem}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{example}[theorem]{Example}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{principle}[theorem]{Principle}
\newenvironment{keywords}{\centerline{\bf\small
Keywords}\begin{quote}\small}{\par\end{quote}\vskip 1ex}

\newtheorem{myexample}[theorem]{Example}
\def\fexample#1#2#3{\vspace{-1ex}\begin{myexample}[#2]\label{#1}\rm #3
\hspace*{\fill} $\diamondsuit\quad$ \end{myexample}\vspace{-2ex} }
\def\paradot#1{\vspace{\paravsp plus 0.5\paravsp minus 0.5\paravsp}\noindent{\bf\boldmath{#1.}}}
\def\paranodot#1{\vspace{\paravsp plus 0.5\paravsp minus 0.5\paravsp}\noindent{\bf\boldmath{#1}}}
\def\req#1{(\ref{#1})}
\def\toinfty#1{\smash{\stackrel{#1\to\infty}{\longrightarrow}}}
\def\iffs{\Leftrightarrow}
\def\eps{\varepsilon}
\def\epstr{\epsilon}                    % for empty string
\def\nq{\hspace{-1em}}
\def\qed{\hspace*{\fill}\rule{1.4ex}{1.4ex}$\quad$\\}
\def\fr#1#2{{\textstyle{#1\over#2}}}
\def\frs#1#2{{^{#1}\!/\!_{#2}}}
\def\SetR{I\!\!R}
\def\SetN{I\!\!N}
\def\SetB{I\!\!B}
\def\SetZ{Z\!\!\!Z}
\def\qmbox#1{{\quad\mbox{#1}\quad}}
\def\e{{\rm e}}                        % natural e
\def\B{\{0,1\}}
\def\E{{\bf E}}                         % Expectation
\def\P{{\rm P}}                         % Probability
\def\lb{{\log_2}}
\def\v{\boldsymbol}
\def\trp{{\!\top\!}}
\def\EE{I\!\!E}
\def\G{\Gamma}
\def\a{\alpha}
\def\d{\delta}
\def\g{\gamma}
\def\s{\sigma}
\def\t{\theta}

\def\b{\beta}
\def\l{\lambda}
\def\r{\rho}

\def\D{{\cal D}}
\def\F{{\cal F}}
\def\M{{\cal M}}
\def\N{{\cal N}}
\def\R{{\cal R}}
\def\X{{\cal X}}
\def\Y{{\cal Y}}

\def\L{\text{\rm Loss}}
\def\h{\hat}

\def\Rank{\text{\rm Rank}}
\def\LR{\text{\rm LR}}
\def\tr{\text{\rm tr}}
\def\Gauss{\text{Gauss}}
\def\KL{\text{\rm KL}}

\begin{document}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                    T i t l e - P a g e                    %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\title{\bf\Large\hrule height5pt \vskip 4mm
The Loss Rank Principle for Model Selection
\vskip 4mm \hrule height2pt}
\author{{\bf Marcus Hutter}\\[3mm]
\normalsize RSISE$\,$@$\,$ANU and SML$\,$@$\,$NICTA \\
\normalsize Canberra, ACT, 0200, Australia \\
\normalsize \texttt{marcus@hutter1.net \ \  www.hutter1.net}
}
\date{February 2006}
\maketitle
\vspace*{-5ex}

\begin{abstract}\noindent
We introduce a new principle for model selection in regression and
classification.
%
Many regression models are controlled by some smoothness or
flexibility or complexity parameter $c$, e.g.\ the number of
neighbors to be averaged over in k nearest neighbor (kNN) regression
or the polynomial degree in regression with polynomials. Let $\h f_D^c$
be the (best) regressor of complexity $c$ on data $D$.
%
A more flexible regressor can fit more data $D'$ well than a more
rigid one. If something (here small loss) is easy to achieve it's
typically worth less. We define the loss rank of $\h f_D^c$ as the
number of other (fictitious) data $D'$ that are fitted better by
$\h f_{D'}^c$ than $D$ is fitted by $\h f_D^c$. We suggest selecting the
model complexity $c$ that has minimal loss rank (LoRP).
%
Unlike most penalized maximum likelihood variants (AIC,BIC,MDL),
LoRP only depends on the regression functions and the loss function.
It works without a stochastic noise model, and is directly
applicable to any non-parametric regressor, like kNN.
%
In this paper we formalize, discuss, and motivate LoRP, study it for
specific regression problems, in particular linear ones, and compare
it to other model selection schemes.
%
\def\contentsname{\centering\normalsize Contents}
{\parskip=-2.7ex\tableofcontents}
\end{abstract}

\vspace*{-2ex}
\begin{keywords}
Model selection,
loss rank principle,
non-parametric regression,
classification
general loss function,
k nearest neighbors.
\end{keywords}

\newpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}\label{secIntro}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%-------------------------------%
\paradot{Regression}
%-------------------------------%
Consider a regression or classification problem in which we want to
determine the functional relationship $y_i\approx f_{true}(x_i)$
from data $D = \{(x_1,y_1),...,(x_n,y_n)\}\in\D$, i.e.\ we seek a
function $f_D$ such that $f_D(x)$ is close to the unknown
$f_{true}(x)$ for all $x$. One may define regressor $f_D$ directly, e.g.\
`average the $y$ values of the $k$ nearest neighbors (kNN) of $x$
in $D$', or select the $f$ from a class of functions $\F$ that has
smallest (training) error on $D$. If the class $\F$ is not too
large, e.g. the polynomials of fixed reasonable degree $d$, this
often works well.

%-------------------------------%
\paradot{Model selection}
%-------------------------------%
What remains is to select the right model complexity $c$, like $k$
or $d$. This selection cannot be based on the training error, since
the more complex the model (large $d$, small $k$) the better the fit
on $D$ (perfect for $d=n$ and $k=1$). This problem is called
overfitting, for which various remedies have been suggested:

We will not discuss empirical test set methods like
cross-validation, but only training set based methods. See e.g.\
\cite{MacKay:92} for a comparison of cross-validation with Bayesian
model selection. Training set based model
selection methods allow using all data $D$ for regression. The
most popular ones can be regarded as penalized versions of Maximum
Likelihood (ML). In addition to the function class $\F$, one has to
specify a sampling model $\P(D|f)$, e.g.\ that the $y_i$ have
independent Gaussian distribution with mean $f(x_i)$. ML chooses $\h
f_D^c=\arg\max_{f\in\F_c}\P(D|f)$, Penalized ML (PML) then chooses
$\h c=\arg\min_c\{-\log\P(D|\h f_D^c)+$Penalty$(c)\}$, where the
penalty depends on the used approach (MDL \cite{Rissanen:78}, BIC
\cite{Schwarz:78}, AIC \cite{Akaike:73}). In particular, modern MDL
\cite{Gruenwald:04} has sound exact foundations and works very well
in practice. All PML variants rely on a proper sampling model (which
may be difficult to establish), ignore (or at least do not tell how
to incorporate) a potentially given loss function, and are typically
limited to (semi)parametric models.

%-------------------------------%
\paradot{Main idea}
%-------------------------------%
The main goal of the paper is to establish a criterion for selecting
the ``best'' model complexity $c$ %
based on regressors $\h f_D^c$ given as a black box without insight into the
origin or inner structure of $\h f_D^c$, %
that does not depend on things often not given (like a stochastic noise model), %
and that exploits what is given (like the loss function).
%
The key observation we exploit is that large classes $\F_c$ or more
flexible regressors $\h f_D^c$ can fit more data $D'\in\D$ well than
more rigid ones, e.g.\ many $D'$ can be fit well with high order
polynomials. We define the {\em loss rank} of $\h f_D^c$ as the number of
other (fictitious) data $D'\in\D$ that are fitted better by
$f_{D'}^c$ than $D$ is fitted by $\h f_D^c$, as measured by some loss
function.
The loss rank is large for regressors fitting $D$ not well {\em and} for
too flexible regressors (in both cases the regressor fits many other
$D'$ better). The loss rank has a minimum for not too flexible
regressors which fit $D$ not too bad. We claim that minimizing
the loss rank is a suitable model selection criterion, since it
trades off the quality of fit with the flexibility of the model.
%
Unlike PML, our new Loss Rank Principle (LoRP) works without a noise
(stochastic sampling) model, and is directly applicable to any
non-parametric regressor, like kNN.

%-------------------------------%
\paradot{Contents}
%-------------------------------%
In Section \ref{secLoRP}, after giving a brief introduction to
regression, we formally state LoRP for model selection.
To make it applicable to real problems, we have to generalize it
to continuous spaces and regularize infinite loss ranks.
%
In Section \ref{secLM} we derive explicit expressions for the loss
rank for the important class of linear regressors, which includes
kNN, polynomial, linear basis function (LBFR), Kernel, and
projective regression.
%
In Section \ref{secBayes} we compare linear LoRP to Bayesian model
selection for linear regression with Gaussian noise and prior, and
%
in Section \ref{secOMS} to PML, in particular MDL, BIC, AIC, and
MacKay's \cite{MacKay:92} and Hastie's et al.\ \cite{Hastie:01}
trace formulas for the effective dimension.
%
In this paper we just scratch at the surface of LoRP. Section
\ref{secMisc} contains further considerations, to be elaborated on
in the future.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{The Loss Rank Principle (LoRP)}\label{secReg}\label{secLoRP}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

After giving a brief introduction to regression, classification,
model selection, overfitting, and some reoccurring examples
(polynomial regression Example \ref{exPR} and kNN Example
\ref{exKNN}), we state our novel Loss Rank Principle for model
selection. We first state it for classification (Principle
\ref{prLRD} for discrete values), and then generalize it for
regression (Principle \ref{prLRC} for continuous values), and exemplify
it on two (over-simplistic) artificial Examples \ref{exSD} and
\ref{exSC}. Thereafter we show how to regularize LoRP for realistic
regression problems.

%-------------------------------%
\paradot{Setup}
%-------------------------------%
We assume data $D = (\v x,\v y) := \{(x_1,y_1),...,(x_n,y_n)\} \in
(\X\times\Y)^n=:\D$ has been observed. We think of the $y$ as having an
approximate functional dependence on $x$, i.e.\ $y_i\approx
f_{true}(x_i)$, where $\approx$ means that the $y_i$ are distorted
by noise or otherwise from the unknown ``true'' values
$f_{true}(x_i)$.

%-------------------------------%
\paradot{Regression and classification}
%-------------------------------%
In regression problems $\Y$ is typically (a subset of) the real
numbers $\SetR$ or some more general measurable space like
$\SetR^m$. In classification, $\Y$ is a finite set or at least
discrete. We impose no restrictions on $\X$. Indeed, $\v x$ will
essentially be fixed and plays only a spectator role, so we will
often notationally suppress dependencies on $\v x$.
%
The goal of regression is to find a function
$f_D\in\F\subset\X\to\Y$ ``close'' to $f_{true}$ based on the past
observations $D$. Or phrased in another way: we are interested in a
regression function $r:\D\to\F$ such that $\h y:=r(x|D)\equiv
r(D)(x)\equiv f_D(x)\approx f_{true}(x)$ for all $x\in\X$.

%-------------------------------%
\paradot{Notation}
%-------------------------------%
We will write $(x,y)$ or $(x_0,y_0)$ for generic data points, use
vector notation $\v x=(x_1,...,x_n)^\trp$ and $\v y=(y_1,...,y_n)^\trp$,
and $D'=(\v x',\v y')$ for generic (fictitious) data of size $n$.

%-------------------------------%
\fexample{exPR}{polynomial regression}{
%-------------------------------%
For $\X=\Y=\SetR$, consider the set $\F_d:=\{f_{\v w}(x)=w_d
x^{d-1}+...w_2 x+w_1 : \v w\in\SetR^d \}$ of polynomials of degree
$d-1$. Fitting the polynomial to data $D$, e.g.\ by least squares
regression, we estimate $\v w$ with $\v{\h w}_D$. The regression
function $\h y=r_d(x|D)=f_{\v{\h w}_D}(x)$ can be written down in
closed form (see Example \ref{exLBFR}).
} % End of Example

%-------------------------------%
\fexample{exKNN}{k nearest neighbors, kNN}{
%-------------------------------%
Let $\Y$ be some vector space like $\SetR$ and $\X$ be a metric
space like $\SetR^m$ with some (e.g.\ Euclidian) metric
$d(\cdot,\cdot)$. kNN estimates $f_{true}(x)$ by averaging the
$y$ values of the $k$ nearest neighbors ${\N}_k(x)$ of $x$ in $D$, i.e.\
$r_k(x|D)={1\over k}\sum_{i\in\N_k(x)} y_i$ %
with $|\N_k(x)|=k$ such that $d(x,x_i)\leq d(x,x_j)$ %
for all $i\in\N_k(x)$ and $j\not\in\N_k(x)$.
} % End of Example

%-------------------------------%
\paradot{Parametric versus non-parametric regression}
%-------------------------------%
Polynomial regression is an example of parametric regression in the
sense that $r_d(D)$ is the optimal function from a family of
functions $\F_d$ indexed by $d<\infty$ real parameters ($\v w$). In
contrast, the kNN regressor $r_k$ is directly given and is not based
on a finite-dimensional family of functions. In general, $r$ may be
given either directly or be the result of an optimization process.

%-------------------------------%
\paradot{Loss function}
%-------------------------------%
The quality of fit to the data is usually measured by a loss function
$\L(\v y,\v{\h y})$, where $\h y_i=\h f_D(x_i)$ is an estimate of $y_i$.
Often the loss is additive: $\L(\v y,\v{\h y})=\sum_{i=1}^n\L(y_i,\h
y_i)$. If the class $\F$ is not too large, good regressors $r$ can
be found by minimizing the loss w.r.t.\ all $f\in\F$. For instance,
$r_d(D)=\arg\min_{f\in\F_d}\sum_{i=1}^n (y_i-f(x_i))^2$ and $\h
y=r_d(x|D)$ in Example
\ref{exPR}.

%-------------------------------%
\paradot{Regression class and loss}
%-------------------------------%
In the following we assume a (typically countable) class of
regressors $\R$ (whatever their origin), e.g.\ the kNN regressors
$\{r_k:k\in\SetN\}$ or the least squares polynomial
regressors $\{r_d:d\in\SetN_0\}$. Note that unlike $f\in\F$,
regressors $r\in\R$ are not functions of $x$ alone but depend
on all observations $D$, in particular on $\v y$.
%
Like for functions $f$, we can compute the loss of each regressor
$r\in\R$:
\beqn
  \L_r(D) \;\equiv\; \L_r(\v y|\v x) \;:=\; \L(\v y,\v{\h y})
  \;=\; \sum_{i=1}^n \L(y_i,r(x_i|\v x,\v y))
\eeqn
where $\h y_i=r(x_i|D)$ in the third expression, and the last
expression holds in case of additive loss.

%-------------------------------%
\paradot{Overfitting}
%-------------------------------%
Unfortunately, minimizing $\L_r$ w.r.t.\ $r$ will typically {\em
not} select the ``best'' overall regressor. This is the well-known
overfitting problem. In case of polynomials, the classes
$\F_d\subset\F_{d+1}$ are nested, hence $\L_{r_d}$ is monotone
decreasing in $d$ with $\L_{r_n}\equiv 0$ perfectly fitting the
data. In case of kNN, $\L_{r_k}$ is more or less an increasing
function in $k$ with perfect regression on $D$ for $k=1$, since no
averaging takes place.
%
In general, $\R$ is often indexed by a ``flexibility'' or smoothness
or complexity parameter, which has to be properly determined. More
flexible $r$ can closer fit the data and hence have smaller loss,
but are not necessarily better, since they have higher variance.
Clearly, too inflexible $r$ also lead to a bad fit (``high bias'').

%-------------------------------%
\paradot{Main goal}
%-------------------------------%
The main goal of the paper is to establish a selection criterion for
the ``best'' regressor $r\in\R$
\begin{itemize}\parskip=0ex\parsep=0ex\itemsep=0ex
\item based on $r$ given as a black box that does not require insight into the
origin or inner structure of $r$,
\item that does not depend on things often not given (like a stochastic noise model),
\item that exploits what is given (like the loss function).
\end{itemize}
%
While for parametric (e.g.\ polynomial) regression, MDL and Bayesian
methods work well (effectively the number of parameters serve as
complexity penalty), their use is seriously limited for
non-parametric black box $r$ like kNN or if a stochastic/coding
model is hard to establish (see Section \ref{secBayes} for a
detailed comparison).

%-------------------------------%
\paradot{Main idea: loss rank}
%-------------------------------%
The key observation we exploit is that a more flexible $r$ can fit
more data $D'\in\D$ well than a more rigid one. For instance, $r_d$
can perfectly fit all $D'$ for $d=n$, all $D'$ that lie on a
parabola for $d=3$, but only linear $D'$ for $d=2$. We consider
discrete $\Y$ i.e.\ classification first, and fix $\v x$. $\v y$ is
the observed data and $\v y'$ are fictitious others.

Instead of minimizing the unsuitable $\L_r(\v y|\v x)$ w.r.t.\ $r$,
we could ask how many $\v y'\in\Y^n$ lead to smaller $\L_r$ than $\v
y$. Many $\v y'$ have small loss for flexible $r$, and so smallness
of $\L_r$ is less significant than if $\v y$ is among very few other
$\v y'$ with small $\L_r$. We claim that the loss rank of $\v y$
among all $\v y'\in\Y^n$ is a suitable measure of fit. We define the
rank of $\v y$ under $r$ as the number of $\v y'\in\Y^n$ with
smaller or equal loss than $\v y$:
\bqa\label{eqRank}
  \Rank_r(\v y|\v x) \;\equiv\; \Rank_r(L)
  &:=& \#\{\v y'\in\Y^n : \L_r(\v y'|\v x)\leq L\},
\\ \nonumber
  \qmbox{where} L &:=& \L_r(\v y|\v x)
\eqa
For this to make sense, we have to assume (and will later assure)
that $\Rank_r(L)<\infty$, i.e.\ there are only
finitely many $\v y'\in\Y^n$ having loss smaller than $L$.
%
In a sense, $\rho=\Rank_r(\v y|\v
x)$ measures how compatible $\v y$ is with $r$; $\v y$ is the
$\rho$th most compatible with $r$.

Since the logarithm is a strictly monotone increasing function, we
can also consider the logarithmic rank $\LR_r(\v y|\v
x):=\log\Rank_r(\v y|\v x)$, which will be more convenient.

%-------------------------------%
\begin{principle}[loss rank principle (LoRP) for classification]\label{prLRD}
%-------------------------------%
For discrete $\Y$, the best classifier/regressor $r:\D\times\X\to\Y$
in some class $\R$ for data $D=(\v x,\v y)$ is the one of smallest
loss rank:
\beq\label{eqLRD}
  r^{best} \;=\; \arg\min_{r\in\R} \LR_r(\v y|\v x)
  \;\equiv\; \arg\min_{r\in\R} \Rank_r(\v y|\v x)
\eeq
where $\Rank_r$ is defined in \req{eqRank}.
\end{principle}

We give a simple example for which we can compute all ranks by hand
to help better grasping how the principle works, but the example is
too simplistic to allow any conclusion on whether the principle is
appropriate.

%-------------------------------%
\fexample{exSD}{simple discrete}{
%-------------------------------%
Consider $\X=\{1,2\}$, $\Y=\{0,1,2\}$, and two points
$D=\{(1,1),(2,2)\}$ lying on the diagonal $x=y$, with polynomial
(zero, constant, linear) least squares regression
$\R=\{r_0,r_1,r_2\}$ (see Ex.\ref{exPR}). $r_0$ is simply 0, $r_1$
the $y$-average, and $r_2$ the line through points $(1,y_1)$ and
$(2,y_2)$. This, together with the quadratic $\L$ for generic $\v
y'$ and observed $\v y=(1,2)$ (and fixed $\v x=(1,2)$), is
summarized in the following table
\beqn
\begin{array}{c|c|c|c}
  d & r_d(x|\v x,\v y') & \L_d(\v y'|\v x) & \L_d(D) \\ \hline
  0 &       0        & y'_1\!\,^2+y'_2\!\,^2 & 5         \\
  1 & \fr12(y'_1+y'_2)  & \fr12(y'_2-y'_1)^2  & \fr12       \\
  2 & (y'_2-y'_1)(x-1)+y'_1 & 0             & 0
\end{array}
\eeqn
From the $\L$ we can easily compute the Rank for all nine $\v
y'\in\{0,1,2\}^2$. Equal rank due to equal loss is
indicated by a $=$ in the table below. Whole equality groups are
actually assigned the rank of their right-most member, e.g.\ for
$d=1$ the ranks of $(y'_1,y'_2)=(0,1),(1,0),(2,1),(1,2)$ are all 7 (and not
4,5,6,7).
\beqn\def\trq{\hspace{4.2ex}}
\begin{array}{c|c|c}
    & \Rank_{r_d}(y'_1y'_2|12) \\
  d & \quad\trq 1 \trq 2 \trq 3 \trq 4 \trq 5 \trq 6 \trq 7 \trq 8 \trq 9 \!\! & \Rank_{r_d}(D)\\ \hline
  0 & y'_1 y'_2 = 00<01=10<11<02=20<21={\bf 12}<22                             & 8             \\
  1 & y'_1 y'_2 = 00=11=22<01=10=21={\bf 12}<02=20                             & 7             \\
  2 & y'_1 y'_2 = 00=01=02=10=11=20=21=22={\bf 12}                             & 9             \\
\end{array}
\eeqn
So LoRP selects $r_1$ as best regressor, since it has minimal rank
on $D$. $r_0$ fits $D$ too badly and $r_2$ is too flexible (perfectly
fits all $D'$).
} % End of Example

%-------------------------------%
\paradot{LoRP for continuous $\Y$}
%-------------------------------%
We now consider the case of continuous or measurable spaces $\Y$,
i.e.\ normal regression problems. We assume $\Y=\SetR$ in the
following exposition, but the idea and resulting principle hold for
more general measurable spaces like $\SetR^m$. We simply reduce the
model selection problem to the discrete case by considering the
discretized space $\Y_\eps=\eps\SetZ$ for small $\eps>0$ and
discretize $\v y\leadsto \v y_\eps\in\eps\SetZ^n$. Then
$\Rank_r^\eps(L):=\#\{\v y'_\eps\in\Y_\eps^n : \L_r(\v y'_\eps|\v
x)\leq L\}$ with $L=\L_r(\v y_\eps|\v x)$ counts the number of
$\eps$-grid points in the set
\beq\label{defV}
  V_r(L) \;:=\; \{\v y'\in\Y^n : \L_r(\v y'|\v x)\leq L\}
\eeq
which we assume (and later assure) to be finite, analogous to the
discrete case. Hence $\Rank_r^\eps(L)\cdot\eps^n$ is an
approximation of the {\em loss volume} $|V_r(L)|$ of set $V_r(L)$, and
typically
$\Rank_r^\eps(L)\cdot\eps^n=|V_r(L)|\cdot(1+O(\eps))\to|V_r(L)|$ for
$\eps\to 0$. Taking the logarithm we get $\LR_r^\eps(\v y|\v
x)=\log\Rank_r^\eps(L)=\log|V_r(L)|-n\log\eps +O(\eps)$. Since
$n\log\eps$ is independent of $r$, we can drop it in comparisons
like \req{eqLRD}. So for $\eps\to 0$ we can define the log-loss
``rank'' simply as the log-volume
\beq\label{defLRC}
  \LR_r(\v y|\v x) \;:=\; \log|V_r(L)|,
  \qmbox{where} L:=\L_r(\v y|\v x)
\eeq

%-------------------------------%
\begin{principle}[loss rank principle for regression]\label{prLRC}
%-------------------------------%
For measurable $\Y$, the best regressor $r:\D\times\X\to\Y$ in some
class $\R$ for data $D=(\v x,\v y)$ is the one of smallest loss
volume:
\beqn\label{eqLRC}
  r^{best} \;=\; \arg\min_{r\in\R} \LR_r(\v y|\v x)
  \;\equiv\; \arg\min_{r\in\R} |V_r(L)|
\eeqn
where $\LR$, $V_r$, and $L$ are defined in \req{defV} and \req{defLRC},
and $|V_r(L)|$ is the volume of $V_r(L)\subseteq\Y^n$.
\end{principle}

For discrete $\Y$ with counting measure we recover the discrete Loss Rank
Principle \ref{prLRD}.

%-------------------------------%
\fexample{exSC}{simple continuous}{
%-------------------------------%
Consider Example \ref{exSD} but with interval $\Y=[0,2]$.
The first table remains unchanged, while the second table becomes
\beqn
\begin{array}{c|c|c|c}
  d & V_d(L)=\{\v y'\in[0,2]^2: ...\} & |V_d(L)| & |V_d(\L_d(D))| \\ \hline
  0 & y'_1\!\,^2+y'_2\!\,^2\leq L  & {{2\sqrt{\max\{L-4,0\}\!}\;+}\atop {L({\pi\over 4}-\cos^{-1}(\min\{ {2\over\sqrt{L}},1\}))}} & \approx 3.6 \\
  1 & \fr12(y'_2-y'_1)^2\leq L      & 4\sqrt{2L}-2L & 3 \\
  2 & 0\leq L                      & 4 & 4
\end{array}
\eeqn
So LoRP again selects $r_1$ as best regressor, since it has smallest loss volume
on $D$.
} % End of Example

%-------------------------------%
\paradot{Infinite rank or volume}
%-------------------------------%
Often the loss rank/volume will be infinite, e.g.\ if we had chosen
$\Y=\SetZ$ in Ex.\ref{exSD} or $\Y=\SetR$ in Ex.\ref{exSC}. We will
encounter such infinities in Section \ref{secLM}. There are various
potential remedies. We could
modify (a) the regressor $r$ or %
(b) the $\L$ to make $\LR_r$ finite, %
(c) the Loss Rank Principle itself, or %
(d) find problem-specific solutions. %
Regressors $r$ with infinite rank might be rejected for
philosophical or pragmatic reasons. We will briefly consider (a) for
linear regression later, but to fiddle around with $r$ in a generic
(blackbox way) seems difficult. We have no good idea how to tinker
with LoRP (c), and also a patched LoRP may be less attractive. For
kNN on a grid we later use remedy (d). While in (decision) theory,
the application's goal determines the loss, in practice the loss is
often more determined by convenience or rules of thumb. So the $\L$
(b) seems the most inviting place to tinker with. A very simple
modification is to add a small penalty term to the loss.
\beq\label{eqLa}
   \L_r(\v y|\v x) \;\leadsto\;
   \L_r^\a(\v y|\v x) := \L_r(\v y|\v x)+\a||\v y||^2,
   \quad \a>0 \mbox{ ``small''}
\eeq
The Euclidian norm $||\v y||^2:=\sum_{i=1}^n y_i^2$ is default, but
other (non)norm regularizes are possible. The regularized
$\LR_r^\a(\v y|\v x)$ based on $\L_r^\a$ is always finite, since
$\{\v y:||\v y||^2\leq L\}$ has finite volume.
%
An alternative penalty $\a\v{\h y}^\trp\v{\h y}$, quadratic in
the regression estimates $\h y_i=r(x_i|\v x,\v y)$ is possible if
$r$ is unbounded in every $\v y\to\infty$ direction.

A scheme trying to determine a single (flexibility) parameter (like $d$
and $k$ in the above examples) would be of no use if it depended
on one (or more) other unknown parameters ($\a$), since varying through
the unknown parameter leads to any (non)desired result.
%
Since LoRP seeks the $r$ of smallest rank, it is natural to also
determine $\a$ by minimizing $\LR_r^\a$ w.r.t.\ $\a$. The good news
is that this leads to meaningful results.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{LoRP for Linear Models}\label{secLM}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

In this section we consider the important class of linear regressors
with quadratic loss function. Since linearity is only assumed in $y$
and the dependence on $x$ can be arbitrary, this class is richer
than it may appear. It includes kNN (Example \ref{exKNN2}), kernel
(Example \ref{exKern}), and many other regressors. For linear
regression and $\Y=\SetR$, the loss rank is the volume of an
$n$-dimensional ellipsoid, which can efficiently be computed in time
$O(n^3)$ (Theorem \ref{thmLRL}). For the special case of projective
regression, e.g.\ linear basis function regression (Example
\ref{exLBFR}), we can even determine the regularization parameter
$\a$ analytically (Theorem \ref{thmLRP}).

%-------------------------------%
\paradot{Linear regression}
%-------------------------------%
We assume $\Y=\SetR$ in this section; generalization to $\SetR^m$ is
straightforward. A linear regressor $r$ can be written in the form
\beq\label{eqmj}
  \h y \;=\; r(x|\v x,\v y) \;=\; \sum_{j=1}^n m_j(x,\v x)y_j
   \quad\forall x\in\X
   \qmbox{and some} m_j:\X\times\X^n\to\SetR
\eeq
Particularly interesting is $r$ for $x=x_1,...,x_n$.
\beq\label{eqM}
  \h y_i \;=\; r(x_i|\v x,\v y) \;=\; \sum_j M_{ij}(\v x)y_j
  \qmbox{with} M:\X^n\to\SetR^{n\times n}
\eeq
where matrix $M_{ij}(\v x)=m_j(x_i,\v x)$. Since LoRP needs $r$ only on the
training data $\v x$, we only need $M$.

%-------------------------------%
\fexample{exKNN2}{kNN ctd.}{
%-------------------------------%
For kNN of Ex.\ref{exKNN} we have
$m_j(x,\v x)={1\over k}$ if $j\in\N_k(x)$ and 0 else, and
$M_{ij}(\v x)={1\over k}$ if $j\in\N_k(x_i)$ and 0 else.
} % End of Example

%-------------------------------%
\fexample{exKern}{kernel regression}{
%-------------------------------%
Kernel regression takes a weighted average over $\v y$,
where the weight of $y_j$ to $y$ is proportional to
the similarity of $x_j$ to $x$, measured by a kernel
$K(x,x_j)$, i.e.\ $m_j(x,\v x)=K(x,x_j)/\sum_{j=1}^n K(x,x_j)$.
For example the Gaussian kernel for $\X=\SetR^m$ is
$K(x,x_j)=\e^{-||x-x_j||_2^2/2\s^2}$.
} % End of Example

%-------------------------------%
\fexample{exLBFR}{linear basis function regression, LBFR}{
%-------------------------------%
Let $\phi_1(x),...,\phi_d(x)$ be a set or vector of ``basis''
functions often called ``features''. We place no restrictions on
$\X$ or $\v\phi:\X\to\SetR^d$. Consider the class of functions
linear in $\v\phi$:
\beqn
  \F_d \;=\; \{ f_{\v w}(x)=\textstyle{\sum_{a=1}^d} w_a\phi_a(x)=\v w^\trp\v\phi(x) : \v w\in\SetR^d \}
\eeqn
For instance, for $\X=\SetR$ and $\phi_a(x)=x^{a-1}$ we would recover
the polynomial regression Example \ref{exPR}.
For quadratic loss function $\L(y_i,\h y_i)=(y_i-\h y_i)^2$ we have
\beqn
  \L_{\v w}(\v y|\v\phi)
  \;:=\; \sum_{i=1}^n(y_i-f_{\v w}(x_i))^2
  \;=\; \v y^\trp\v y - 2\v y^\trp\v\Phi\v w + \v w^\trp B \v w
\eeqn
where matrix $\v\Phi$ is defined by $\v\Phi_{ia}=\phi_a(x_i)$ and
$B$ is a symmetric matrix with
$B_{ab}=\sum_{i=1}^n\phi_a(x_i)\phi_b(x_i)=[\v\Phi^\trp\v\Phi]_{ab}$.
The loss is quadratic in $\v w$ with minimum at $\v
w=B^{-1}\v\Phi^\trp\v y$. So the least squares regressor is $\h y=\v
y^\trp\v\Phi B^{-1}\v\phi(x)$, hence $m_j(x,\v x)=(\v\Phi
B^{-1}\v\phi(x))_j$ and $M(\v x)=\v\Phi B^{-1}\v\Phi^\trp$.
} % End of Example

Consider now a general linear regressor $M$ with quadratic loss
and quadratic penalty
\bqa\nonumber
  \L_M^\a(\v y|\v x) &=& \sum_{i=1}^n
  \left(y_i-\textstyle{\sum_{j=1}^n} M_{ij}y_j\right)^2+\a||\v y||^2
  \;=\; \v y^\trp S_\a\v y,
\\ \label{defSa}
  \qmbox{where\footnotemark} S_\a &=& (1\!\!1-M)^\trp(1\!\!1-M)+\a 1\!\!1
\eqa%
\footnotetext{The mentioned alternative penalty $\a||\v{\h
y}||^2$ would lead to
$S_\a = (1\!\!1-M)^\trp(1\!\!1-M)+\a M^\trp M$.
For LBFR, penalty $\a||\v{\h w}||^2$ is popular (ridge regression).
Apart from being limited to parametric regression, it
has the disadvantage of not being reparametrization invariant.
For instance, scaling $\phi_a(x)\leadsto\g_a\phi_a(x)$ doesn't
change the class $\F_d$, but changes the ridge regressor.}
%
($1\!\!1$ is the identity matrix). $S_\a$ is a symmetric matrix. For
$\a>0$ it is positive definite and for $\a=0$ positive semidefinite.
If $\l_1,...\l_n\geq 0$ are the eigenvalues of $S_0$, then
$\l_i+\a$ are the eigenvalues of $S_\a$. $V(L)=\{\v y'\in\SetR^n : \v
y'\!\,^\trp S_\a\v y'\leq L\}$ is an ellipsoid with the eigenvectors of
$S_\a$ being the main axes and $\sqrt{L/(\l_i+\a)}$ being their length.
Hence the volume is
\beqn
  |V(L)| \;=\; v_n\prod_{i=1}^n \sqrt{L\over \l_i+\a}
  \;=\; {v_n L^{n/2}\over\sqrt{\det S_\a}}
\eeqn
where $v_n=\pi^{n/2}/{n\over 2}!$ is the volume of the
$n$-dimensional unit sphere, $z!:=\G(z+1)$, and $\det$ is the
determinant. Taking the logarithm we get
\beq\label{eqLRL1}
  \LR_M^\a(\v y|\v x)
  \;=\; \log|V(\L_M^\a(\v y|\v x))|
  \;=\; \fr n2\log(\v y^\trp S_\a\v y)-\fr12\log\det S_\a +\log v_n
\eeq
Consider now a {\em class} of linear regressors $\R=\{M\}$,
e.g.\ the kNN regressors $\{M_k:k\in\SetN\}$ or
the $d$-dimensional linear basis function regressors $\{M_d:d\in\SetN_0\}$.

%-------------------------------%
\begin{theorem}[LoRP for linear regression]\label{thmLRL}
%-------------------------------%
For $\Y=\SetR$, the best linear regressor $M:\X^n\to\SetR^{n\times n}$
in some class $\M$ for data
$D=(\v x,\v y)$ is
\beq\label{eqLRL}
  M^{best} \;=\; \mathop{\arg\min}_{M\in\M,\a\geq 0}
  \{ \fr n2\log(\v y^\trp S_\a\v y)-\fr12\log\det S_\a \}
  \;=\; \mathop{\arg\min}_{M\in\M\;\a\geq 0}
  \Big\{ {\v y^\trp S_\a\v y\over(\det S_\a)^{1/n}} \Big\}
\eeq
where $S_\a$ is defined in \req{defSa}.
\end{theorem}

Since $v_n$ is independent of $\a$ and $M$ it was possible to drop
$v_n$. The last expression shows that linear LoRP minimizes the $\L$
times the geometric average of the squared axes lengths of ellipsoid
$V(1)$. Note that $M^{best}$ depends on $\v y$ unlike the $M\in\M$.

%-------------------------------%
\paradot{Nullspace of $S_0$}
%-------------------------------%
If $M$ has an eigenvalue 1, then $S_0=(1\!\!1-M)^\trp(1\!\!1-M)$ has a zero
eigenvalue and $\a>0$ is necessary, since $\det S_0=0$. Actually
this is true for most practical $M$. Nearly all linear regressors
are invariant under a constant shift of $\v y$, i.e.\
$r(y_i+c|D)=r(y_i|D)+c$, which implies that $M$ has eigenvector
$(1,...,1)^\trp$ with eigenvalue 1. This can easily be checked for kNN
(Ex.\ref{exKNN}), Kernel (Ex.\ref{exKern}), and LBFR
(Ex.\ref{exLBFR}). Such a generic 1-eigenvector effecting all $M\in\M$
could easily and maybe should be filtered out by considering only
the orthogonal space or dropping these $\l_i=0$ when computing $\det
S_0$. The 1-eigenvectors that depend on $M$ are the ones where we
really need a regularizer $\a>0$ for. For instance, $M_d$ in LBFR has $d$
eigenvalues 1, and $M_{\text{kNN}}$ has as many eigenvalues 1 as there are
disjoint components in the graph determined by the edges $M_{ij}>0$
%
In general we need to find the optimal $\a$ numerically.
If $M$ is a projection we can find $\a_{min}$ analytically.

%-------------------------------%
\paradot{Projective regression}
%-------------------------------%
Consider a projection matrix $M=P=P^2$ with $d=\tr P$ eigenvalues 1,
and $n-d$ zero eigenvalues. For instance, $M=\v\Phi B^{-1}\v\Phi^\trp$
of LBFR Ex.\ref{exLBFR} is such a matrix, since $M\v\Phi=\v\Phi$ and
$M\v\Psi=0$ for $\v\Psi$ such that $\v\Phi^\trp\v\Psi=0$. This implies
that $S_\a$ has $d$ eigenvalues $\a$ and $n-d$ eigenvalues $1+\a$.
Hence
\bqa\nonumber
  \det S_\a &=& \a^d(1+\a)^{n-d}, \qmbox{where}
  S_\a = S_0+\a 1\!\!1 = 1\!\!1-P+\a 1\!\!1
\\ \nonumber
  \v y^\trp S_\a\v y &=& (\r+\a)\v y^\trp\v y, \qmbox{where}
  \r:={\v y^\trp S_0\v y\over\v y^\trp\v y} = 1-{\v y^\trp P\v y\over\v y^\trp\v y}
\\ \label{eqLRPa}
  \Rightarrow\quad \LR_P^\a &=& \fr n2\log\v y^\trp\v y + \fr n2\log(\r+\a)-
  {\textstyle{d\over 2}}\log\a - {\textstyle{n-d\over 2}}\log(1+\a)
\eqa
The first term is independent of $\a$. Consider $1-\r>{d\over n}$,
the reasonable region in practice. Solving $\partial
\LR_P^\a/\partial\a = 0$ w.r.t.\ $\a$ we get a minimum at
$\a=\a_{min}:={\r d\over(1-\r)n-d}$. After some algebra we get
\beq\label{eqLRPamin}\textstyle
  \LR_P^{\a_{min}} \;=\; \fr n2\log\v y^\trp\v y - \fr n2 \KL({d\over n}||1-\r),
  \qmbox{where} \KL(p||q) \;=\; p\log{p\over q}+(1-p)\log{1-p\over 1-q}
\eeq
is the relative entropy or Kullback-Leibler divergence.
%
Minimizing $\LR_P^{\a_{min}}$ w.r.t.\ $M$ is equivalent to
maximizing $\KL({d\over n}||1-\r)$. This is an unusual task,
since one mostly encounters $D$ minimizations.
%
For fixed $d$, $\LR_P^{\a_{min}}$ is monotone increasing in $\r$.
Since $\L_P^\a\propto\r+\a$, LoRP suggests to minimize $\L$ for fixed
model dimension $d$.
%
For fixed $\r$, $\LR_P^{\a_{min}}$ is monotone increasing in $d$,
i.e.\ LoRP suggests to minimize model dimension $d$ for fixed $\L$.
Normally there is a tradeoff between minimizing $d$ and $\r$,
and LoRP suggests that the optimal choice is the one that maximizes $\KL$.

%-------------------------------%
\begin{theorem}[LoRP for projective regression]\label{thmLRP}
%-------------------------------%
The best projective regressor $P:\X^n\to\SetR^{n\times n}$
with $P=P^2$ in some projective class $\cal P$ for data
$D=(\v x,\v y)$ is
\beqn\label{eqLRP}
  P^{best} \;=\; \arg\max_{P\in\cal P} \;\textstyle \KL({\tr P(x)\over n}||{\v y^\trp P(x)\v y\over\v y^\trp\v y}),
  \qmbox{provided} {\tr P\over n} < {\v y^\trp P\v y\over\v y^\trp\v y}
\eeqn
\end{theorem}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Comparison to Gaussian Bayesian Linear Regression}\label{secBayes}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

We now consider linear basis function regression (LBFR) from a
Bayesian perspective with Gaussian noise and prior, and compare it
to LoRP. In addition to the noise model as in PML, one also has to
specify a prior. Bayesian model selection (BMS) proceeds by
selecting the model that has largest evidence. In the special case
of LBFR with Gaussian noise and prior and an ML-II estimate for the
noise variance, the expression for the evidence has a similar
structure as the expression of the loss rank.

%-------------------------------%
\paradot{Gaussian Bayesian LBFR / MAP}
%-------------------------------%
Recall from Sec.\ref{secLM} Ex.\ref{exLBFR} that $\F_d$ is the class
of functions $f_{\v w}(x)=\v w^\trp\v\phi(x)$ ($\v w\in\SetR^d$) that
are linear in feature vector $\v\phi$. Let
\beq\label{eqGauss}
  \Gauss_N(\v z|\v\mu,\Sigma) \;:=\;
  {\exp(-\fr12(\v z-\v\mu)^\trp\Sigma^{-1}(\v z-\v\mu))
   \over(2\pi)^{N/2}\sqrt{\det\Sigma}}
\eeq
denote a general $N$-dimensional Gaussian distribution with mean
$\v\mu$ and covariance matrix $\Sigma$.
%
We assume that observations $y$ are perturbed from $f_{\v w}(x)$ by
independent additive Gaussian noise with variance $\b^{-1}$ and zero
mean, i.e.\ the likelihood of $\v y$ under model $\v w$ is
$\P(\v y|\v w)=\Gauss_n(\v y|\v\Phi\v w,\b^{-1}1\!\!1)$,
where $\v\Phi_{ia}=\v\phi_a(x_i)$.
A Bayesian assumes a prior (before seeing $\v y$) distribution on
$\v w$. We assume a centered Gaussian with covariance matrix $(\a
C)^{-1}$, i.e.\
$\P(\v w)=\Gauss_d(\v w|\v 0,\a^{-1}C^{-1})$.
From the prior and the likelihood one can compute the evidence and the posterior
\bqa\label{eqEv}
  \mbox{Evidence:}\qquad\quad\; \P(\v y) &=& \int\P(\v y|\v w)\P(\v w)d\v w
                                      \;=\; \Gauss_n(\v y|\v 0,\b^{-1}S^{-1})
\\ \nonumber
  \mbox{Posterior:}\qquad \P(\v w|\v y) &=& \P(\v y|\v w)\P(\v w)/P(\v y)
                                       \;=\; \Gauss_d(\v w|\v{\h w},A^{-1})
\eqa
\beq\label{eqBAMS}
  B:=\v\Phi^\trp\v\Phi, \quad
  A:=\a C+\b B, \quad
  M:=\b\v\Phi A^{-1}\v\Phi^\trp, \quad
  S:=1\!\!1-M, \quad
\eeq
\beqn
  \v{\h w}:=\b A^{-1}\v\Phi^\trp\v y, \quad
  \v{\h y}:=\v\Phi\v{\h w}=M\v y
\eeqn
A standard Bayesian point estimate for $\v w$ for fixed $d$ is the
one that maximizes the posterior (MAP) (which in the Gaussian case
coincides with the mean)
$\v{\h w}=\arg\max_{\v w}\P(\v w|\v y)=\b A^{-1}\v\Phi^\trp\v y$.
For $\a\to 0$, MAP reduces to Maximum Likelihood (ML), which
in the Gaussian case coincides with the least squares regression of
Ex.\ref{exLBFR}. For $\a>0$, the regression matrix $M$ is not a
projection anymore.

%-------------------------------%
\paradot{Bayesian model selection}
%-------------------------------%
Consider now a family of models $\{\F_1,\F_2,...\}$. Here the $\F_d$
are the linear regressors with $d$ basis functions, but in general
they could be completely different model classes. All quantities in
the previous paragraph implicitly depend on the choice of $\F$,
which we now explicate with an index. In particular, the evidence
for model class $\F$ is $\P_\F(\v y)$.
%
Bayesian Model Selection (BMS) chooses the model class (here $d$) $\F$ of
highest evidence:
\beqn
  \F^{\text{BMS}} \;=\;\arg\max_\F\P_\F(\v y)
\eeqn
Once the model class $\F^{\text{BMS}}$ is determined, the MAP (or other)
regression function $f_{{\v w}_{\F^{\text{BMS}}}}$ or $M_{\F^{\text{BMS}}}$ are
chosen. The data variance $\b^{-1}$ may be known or estimated
from the data, $C$ is often chosen $1\!\!1$, and $\a$ has to be chosen
somehow. Note that while $\a\to 0$ leads to a reasonable
MAP=ML regressor for fixed $d$, this limit cannot be used for BMS.

%-------------------------------%
\paradot{Comparison to LoRP}
%-------------------------------%
Inserting \req{eqGauss} into \req{eqEv} and taking the
logarithm we see that BMS minimizes
\beq\label{eqLLG1}
  -\log\P_\F(\v y) \;=\; \fr\b2\v y^\trp S\v y - \fr12\log\det S - \fr n2\log\fr\b{2\pi}
\eeq
w.r.t.\ $\F$. Let us estimate $\b$ by ML: We assume a broad prior
$\a\ll\b$ so that $\b{\partial S\over\partial\b}=O({\a\over\b})$ can
be neglected. Then ${\partial\log\P_\F(\v y)\over\partial\b} =
\fr12\v y^\trp S\v y-{n\over 2\b}+O({\a\over\b}n)=0$ $\Leftrightarrow$
$\b\approx \h\b:=n/(\v y^\trp S\v y)$. Inserting $\h\b$
into \req{eqLLG1} we get
\beq\label{eqLLG2}
  -\log\P_\F(\v y) \;=\;
  \fr n2\log\v y^\trp S\v y - \fr12\log\det S - \fr n2\log\fr{n}{2\pi\e}
\eeq
Taking an improper prior $\P(\b)\propto\b^{-1}$ and integrating out
$\b$ leads for small $\a$ to a similar result. The last term in
\req{eqLLG2} is a constant independent of $\F$ and can be ignored.
The first two terms have the same structure as in linear LoRP
\req{eqLRL}, but the matrix $S$ is different.
%
In both cases, $\a$ act as regularizers, so we may minimize over
$\a$ in BMS like in LoRP. For $\a=0$ (which neither makes sense in
BMS nor in LoRP), $M$ in BMS coincides with $M$ of Ex.\ref{exLBFR},
but still the $S_0$ in LoRP is the square of the $S$ in BMS. For $\a>0$, $M$ of
BMS may be regarded as a regularized regressor as suggested in
Sec.\ref{secLoRP} (a), rather than a regularized loss function (b) used
in LoRP. Note also that BMS is limited to (semi)parametric regression,
i.e.\ does not cover the non-parametric kNN Ex.\ref{exKNN} and Kernel
Ex.\ref{exKern}, unlike LoRP.

Since $B$ only depends on $\v x$ (and not on $\v y$), and all $\P$
are implicitly conditioned on $\v x$, one could choose $C=B$. In
this case, $M=\g\v\Phi B^{-1}\v\Phi^\trp$, with $\g={\b\over\a+\b}<1$
for $\a>0$, is a simple multiplicative regularization of projection
$\v\Phi B^{-1}\v\Phi^\trp$, and \req{eqLLG2} coincides with
\req{eqLRPa} for suitable $\a$, apart from an irrelevant additive
constant, hence minimizing \req{eqLLG2} over $\a$
also leads to \req{eqLRPamin}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Comparison to other Model Selection Schemes}\label{secOMS}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

In this section we give a brief introduction to Penalized Maximum
Likelihood (PML) for (semi)parametric regression, and its major
instantiations, the Akaike and the Bayesian Information Criterion
(AIC and BIC), and the Minimum Description Length (MDL) principle,
whose penalty terms are all proportional to the number of parameters
$d$. The {\em effective} number of parameters is often much smaller than
$d$, e.g.\ if there are soft constraints like in ridge regression. We
compare MacKay's \cite{MacKay:92} trace formula for Gaussian
Bayesian LBFR and Hastie's et al.\ \cite{Hastie:01} trace formula
for general linear regression with LoRP.

%-------------------------------%
\paradot{Penalized ML (AIC, BIC, MDL)}
%-------------------------------%
Consider a $d$-dimensional stochastic model class like the Gaussian
Bayesian linear regression example of Section \ref{secBayes}. Let
$\P_d(\v y|\v w)$ be the data likelihood under $d$-dimensional model
$\v w\in\SetR^d$. The maximum likelihood (ML) estimator for fixed
$d$ is
\beqn
  \v{\h w} \;=\; \arg\max_{\v w}\P_d(\v y|\v w)
  \;=\; \arg\min_{\v w}\{-\log \P_d(\v y|\v w)\}
\eeqn
Since $-\log\P_d(\v y|\v w)$ decreases with $d$, we
cannot find the model dimension by simply minimizing over $d$
(overfitting). Penalized ML adds a complexity term to get
reasonable results
\beqn
  \h d \;=\; \arg\min_d\{-\log \P_d(\v y|\v{\h w}) + \mbox{Penalty}(d) \}
\eeqn
The penalty introduces a tradeoff between the first and second
term with a minimum at $\h d<\infty$. Various penalties have been
suggested: The Akaike Information Criterion (AIC) \cite{Akaike:73}
uses $d$, the Bayesian Information Criterion (BIC) \cite{Schwarz:78}
and the (crude) Minimum Description Length (MDL) principle use $\fr
d2\log n$ \cite{Rissanen:78,Gruenwald:04} for Penalty$(d)$.
%
There are at least {\em three important conceptual differences} to LoRP:
\begin{itemize}\parskip=0ex\parsep=0ex\itemsep=0ex
\item In order to apply PML one needs to specify not only a class
of regression functions, but a full probabilistic model $\P_d(\v y|\v w)$,
\item PML ignores or at least does not tell how to incorporate
a potentially given loss-function,
\item PML (AIC,BIC,MDL) is mostly limited to (semi)parametric
models (with $d$ ``true'' parameters).
\end{itemize}

We discuss two approaches to the last item in the remainder of this
section: AIC, BIC, and MDL are not directly applicable %
(a) for non-parametric models like kNN or Kernel regression, or %
(b) if $d$ does not reflect the ``true'' complexity of the model.
For instance, ridge regression can work even for $d$ larger than
$n$, because a penalty pulls most parameters towards (but not
exactly to) zero. %
MacKay \cite{MacKay:92} suggests an expression for the effective
number of parameters $d_{e\!f\!f}$ as a substitute for $d$ in case
of (b), and Hastie et al.\ \cite{Hastie:01} more generally also for
(a).

%-------------------------------%
\paradot{The trace penalty for parametric Gaussian LBFR}
%-------------------------------%
We continue with the Gaussian Bayesian linear regression example
(see Section \ref{secBayes} for details and notation). Performing
the integration in \req{eqEv}, MacKay \cite[Eq.(21)]{MacKay:92}
derives the following expression for the Bayesian evidence for $C=1\!\!1$
\bqa\label{eqLLMK}
  -\log\P(\v y) &=&
  (\a \h E_W+\b \h E_D) + (\fr12\log\det A -\fr d2\log\a) - \fr n2\log\fr\b{2\pi}
\\ \nonumber
  \h E_D &=& \fr12||\v\Phi\v{\h w}-\v y||_2^2, \quad
  \h E_W =\fr12||\v{\h w}||_2^2
\eqa
(the first bracket in \req{eqLLMK} equals $\fr\b2\v y^\trp S\v y$ and
the second equals $-\fr12\log\det S$, cf.\ \req{eqLLG1}).
Minimizing \req{eqLLMK} w.r.t.\ $\a$ leads
to the following relation:
\beqn
  0 \;=\; \textstyle{-\partial\log\P(\v y)\over\partial\a} \;=\;
  \h E_W +\fr12\tr A^{-1}-\fr d{2\a} \qquad
 ({\partial\over\partial\a}\log\det A=\tr A^{-1})
\eeqn
He argues that
$\a||\v{\h w}||_2^2$ corresponds to the effective number of
parameters, hence
\beq\label{eqdeffMcK}
  d^{\text{McK}}_{e\!f\!f} \;:=\; \a||\v{\h w}||_2^2
  \;=\; 2\a \h E_W \;=\; d-\a\tr A^{-1}
\eeq

%-------------------------------%
\paradot{The trace penalty for general linear models}
%-------------------------------%
We now return to general linear regression $\v{\h y}=M(\v x)\v y$
\req{eqM}. LBFR is a special case of a projection matrix $M=M^2$
with rank $d=\tr M$ being the number of basis functions. $M$ leaves
$d$ directions untouched and projects all other $n-d$ directions to
zero. For general $M$, Hastie et al. \cite[Sec.5.4.1]{Hastie:01}
argue to regard a direction that is only somewhat shrunken, say by a
factor of $0<\b<1$, as a fractional parameter ($\b$ degrees of
freedom). If $\b_1,...,\b_n$ are the shrinkages = eigenvalues
of $M$, the effective number of parameters could be defined as
\cite[Sec.7.6]{Hastie:01}
\beqn
  d^{\text{HTF}}_{e\!f\!f} \;:=\; \sum_{i=1}^n \b_i \;=\; \tr M
\eeqn
which generalizes the relation $d=\tr M$ beyond projections.
For MacKay's $M$ \req{eqBAMS}, $\tr M=d-\tr A^{-1}$,
i.e.\ $d^{\text{HTF}}_{e\!f\!f}$ is consistent with and generalizes
$d^{\text{McK}}_{e\!f\!f}$.

%-------------------------------%
\paradot{Problems}
%-------------------------------%
Though nicely motivated, the trace formula is not without problems.
First, since for projections, $M=M^2$, one could equally well have
argued for $d^{\text{HTF}}_{e\!f\!f}=\tr M^2$. Second, for kNN we have $\tr
M=\fr nk$ (since $M$ is $\fr 1k$ on the diagonal), which does not
look unreasonable. Consider now kNN' where we average over the $k$
nearest neighbors {\em excluding} the closest neighbor. For
sufficiently smooth functions, kNN' for suitable $k$ is still a
reasonable regressor, but $\tr M=0$ (since $M$ is zero on the
diagonal). So $d^{\text{HTF}}_{e\!f\!f}=0$ for kNN', which makes no sense
and would lead one to always select the $k=1$ model.

%-------------------------------%
\paradot{Relation to LoRP}
%-------------------------------%
In the case of kNN', $\tr M^2$ would be a better estimate for the
effective dimension. In linear LoRP, $-\log\det S_\a$ serves as
complexity penalty. Ignoring the nullspace of
$S_0=(1\!\!1-M)^\trp(1\!\!1-M)$ \req{defSa}, we can Taylor expand
$-\fr12\log\det S_0$ in $M$
\beqn
  -\fr12\log\det S_0 \;=\; -\tr\log(1\!\!1\!-\!M)
  \;=\; \sum_{s=1}^\infty \fr1s\tr(M^s)
  \;=\; \tr M + \fr12\tr M^2 + ...
\eeqn
For BMS \req{eqLLG2} with $S=1\!\!1-M$ \req{eqBAMS} we get half of
this value. So the trace penalty may be regarded as a leading order
approximation to LoRP. The higher order terms prevent
peculiarities like in kNN'.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Outlook}\label{secMisc}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

So far we have only scratched at the surface of the Loss Rank
Principle. LoRP seems to be a promising principle with a lot of
potential, leading to a rich field. In the following we briefly
summarize miscellaneous considerations, which may be elaborated on
in the future: Experiments, Monte Carlo estimates for non-linear
LoRP, numerical approximation of $\det S_\a$, LoRP for
classification, self-consistent regression, explicit expressions for
kNN on a grid, loss function selection, and others.

%-------------------------------%
\paradot{Experiments}
%-------------------------------%
Preliminary experiments on selecting $k$ in kNN regression confirm
that LoRP selects a ``good'' $k$. (Even on artificial data we cannot
determine whether the ``right'' $k$ is selected, since kNN is not a
generative model). LoRP for LBFR seems to be consistent with rapid
convergence.

%-------------------------------%
\paradot{Monte Carlo estimates for non-linear LoRP}
%-------------------------------%
For non-linear regression we did not present an efficient algorithm
for the loss rank/volume $\LR_r(\v y|\v x)$. The high-dimensional
volume $|V_r(L)|$ \req{defV} may be computed by Monte Carlo
algorithms. Normally $V_r(L)$ constitutes a small part of $\Y^n$, and
uniform sampling over $\Y^n$ is not feasible. Instead one should
consider two competing regressors $r$ and $r'$ and compute $|V\cap
V'|/|V|$ and $|V\cap V'|/|V'|$ by uniformly sampling from $V$ and
$V'$ respectively e.g.\ with a Metropolis-type algorithm. Taking the
ratio we get $|V'|/|V|$ and hence the loss rank difference
$\LR_r-\LR_{r'}$, which is sufficient for LoRP. The usual tricks and
problems with sampling apply here too.

%-------------------------------%
\paradot{Numerical approximation of $\det S_\a$}
%-------------------------------%
Even for linear regression, a Monte Carlo algorithm may be faster
than the naive $O(n^3)$ algorithm \cite{Bai:96}. Often $M$ is a very
sparse matrix (like in kNN) or can be well approximated by a sparse
matrix (like for Kernel regression), which allows to approximate
$\det S_\a$, sometimes in linear time \cite{Reusken:02}.

%-------------------------------%
\paradot{LoRP for classification}
%-------------------------------%
A classification problem is or can be regarded as a regression problem
in which $\Y$ is finite. This implies that we need to compute (count)
$\LR_r$ for non-linear $r$ somehow, e.g.\ approximately by Monte Carlo.

%-------------------------------%
\paradot{Self-consistent regression}
%-------------------------------%
So far we have considered only ``on-data'' regression. LoRP only
depends on the regressor $r$ on data $D$ and not on
$x\not\in\{x_1,...,x_n\}$. One can construct canonical regressors
for off-data $x$ from regressors given only on-data in the following
way:
%
We add a virtual data point $(x,y)$ to $D$, where $x$ is the
off-data point of interest. If we knew $y$ we could estimate $\h
y=r(x|\{(x,y)\}\cup D)$, but we don't know $y$. But if we require
consistency, namely that $\h y=y$, we get a canonical
estimate for $\h y$.
%
First, this bootstrap may ease the specification of the regression
models, second, it is a canonical way for interpolation (LoRP
can't distinguish between $r$ that are identical on $D$), and third,
many standard regressors (kNN, Kernel, LBFR) are self-consistent in
the sense that they are canonical.

%-------------------------------%
\paradot{Explicit expressions for kNN on a grid}
%-------------------------------%
In order to get more insight into LoRP, a case that allows an
analytic solution is useful. For k nearest neighbors classification
with $\v x$ lying on a hypercube of the regular grid $\X=\SetZ^d$
one can derive explicit expressions for the loss rank as a function
of $k$, $n$, and $d$. For $n\gg k\gg 3^d$, the penalty
$-\fr12\log\det S$ is proportional to $\tr M$ with proportionality
constant decreasing from about 3.2 for $d=1$ to 1.5 for
$d\to\infty$.

%-------------------------------%
\paradot{LoRP for hybrid model classes}
%-------------------------------%
LoRP is not restricted to model classes indexed by a
single integral ``complexity'' parameter, but may be applied more
generally to selecting among some (typically discrete) class of
models/regressors. For instance, the class could contain kNN {\em
and} polynomial regressors, and LoRP selects the complexity {\em
and} type of regressor (non-parametric kNN versus parametric
polynomials).

%-------------------------------%
\paradot{General additive loss}
%-------------------------------%
Linear LoRP $\v{\h y}=M(\v x)\v y$ of Section \ref{secLM} can easily
be generalized from quadratic to $\rho$-norm $\L_M(\v y|\v
x)=||\v y-\v{\h y}||_\rho^p$ (any $p$). For $\a=0$, $\v y^\trp S_0\v y$ in
\req{eqLRL1} becomes $||\v y-\v{\h y}||_\rho^2$ and $v_\rho$ the
volume of the unit $d$-dimensional $\rho$-norm ``ball''. Useful
expressions for general additive $\L_N=\sum_i h(y_i-\h y_i)$
can also be derived. Regularization may be performed by
$M\leadsto\gamma M$ with optimization over $\gamma<1$.

%-------------------------------%
\paradot{Loss-function selection}
%-------------------------------%
In principle, the loss function should be part of the problem
specification, since it characterizes the ultimate goal.
In reality, though, having to specify the loss function can be a
nuisance.
%
We could interpret the regularized loss \req{eqLa} as a class of
loss functions parameterized by $\a$, and $\arg\min_\a\LR_r^a$ as a
loss function optimization or selection. This suggests to choose in
general the loss function that has minimal loss rank. This leads to
sensible results if the considered class of loss functions is not
too large (e.g.\ all $\rho$-norm losses in the previous paragraph).
So LoRP can be used not only for model selection, but also for
loss function selection.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%         Bibliography        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\addcontentsline{toc}{section}{\refname}
\begin{small}
\begin{thebibliography}{BFG96}

\bibitem[Aka73]{Akaike:73}
H.~Akaike.
\newblock Information theory and an extension of the maximum likelihood
  principle.
\newblock In {\em Proc. 2nd International Symposium on Information Theory},
  pages 267--281, Budapest, Hungary, 1973. Akademiai Kaid\'o.

\bibitem[BFG96]{Bai:96}
Z.~Bai, M.~Fahey, and G.~Golub.
\newblock Some large-scale matrix computation problems.
\newblock {\em Journal of Computational and Applied Mathematics},
  74(1--2):71--89, 1996.

\bibitem[Gr{\"u}04]{Gruenwald:04}
P.~D. Gr{\"u}nwald.
\newblock Tutorial on minimum description length.
\newblock In {\em Minimum Description Length: recent advances in theory and
  practice}, page Chapters 1 and 2. MIT Press, 2004.
\newblock http://www.cwi.nl/$_{^\sim}$pdg/ftp/mdlintro.pdf.

\bibitem[HTF01]{Hastie:01}
T.~Hastie, R.~Tibshirani, and J.~H. Friedman.
\newblock {\em The Elements of Statistical Learning}.
\newblock Springer, 2001.

\bibitem[Mac92]{MacKay:92}
D.~J.~C. MacKay.
\newblock Bayesian interpolation.
\newblock {\em Neural Computation}, 4(3):415--447, 1992.

\bibitem[Reu02]{Reusken:02}
A.~Reusken.
\newblock Approximation of the determinant of large sparse symmetric positive
  definite matrices.
\newblock {\em SIAM Journal on Matrix Analysis and Applications},
  23(3):799--818, 2002.

\bibitem[Ris78]{Rissanen:78}
J.~J. Rissanen.
\newblock Modeling by shortest data description.
\newblock {\em Automatica}, 14(5):465--471, 1978.

\bibitem[Sch78]{Schwarz:78}
G.~Schwarz.
\newblock Estimating the dimension of a model.
\newblock {\em Annals of Statistics}, 6(2):461--464, 1978.

\end{thebibliography}
\end{small}

\end{document}

%-----------------------End-of-LoRP.tex-----------------------%
