wtls_manual.tex

加权总体最小二乘matlab工具箱

\documentclass[10pt]{article}
%-------------------------------------------------------------
\usepackage{graphicx}
\usepackage{psfrag}
\usepackage{amsmath}
\usepackage{amsthm}
\usepackage{amsfonts}
\usepackage{mathptmx}
\usepackage{moreverb}
\usepackage[all]{xy}
\usepackage{rotating}
%--------------------------------------------------------------
\newcommand{\la}[1]{\label{#1}}
\newcommand{\re}[1]{(\ref{#1})}
\newcommand{\bmx}{\begin{bmatrix}}
\newcommand{\emx}{\end{bmatrix}}
\newcommand{\bsm}{\left[\begin{smallmatrix}}
\newcommand{\esm}{\end{smallmatrix}\right]}
\newcommand{\R}{\ensuremath{\mathbb{R}}} % real numbers
\newcommand{\N}{\ensuremath{\mathbb{N}}} % natural numbers
\newcommand{\C}{\ensuremath{\mathbb{C}}} % complex numbers
\newcommand{\Z}{\ensuremath{\mathbb{Z}}} % integer numbers
% --- Basic linear algebra definitions
\newcommand{\rank}{\operatorname{rank}}
\newcommand{\tr}{\operatorname{trace}}
\newcommand{\diag}{\operatorname{diag}}
\newcommand{\blkdiag}{\operatorname{blk\,diag}}
\newcommand{\abs}[1]{|#1|}
\newcommand{\norm}[1]{\|#1\|}
\newcommand{\mtov}{\operatorname{vec}}
\newcommand{\bo}{\text{\bf 1}} % vector of ones
\newcommand{\F}{\text{\rm F}}  % Frobenius norm
\newcommand{\col}{\operatorname{col}} % column vector
\newcommand{\coldim}{\operatorname{col\,dim}}
\newcommand{\rowdim}{\operatorname{row\,dim}}
\newcommand{\image}{\operatorname{image}}
\newcommand{\colspan}{\operatorname{col\,span}} 
\newcommand{\rowspan}{\operatorname{row\,span}}
\newcommand{\degree}{\operatorname{degree}}
\newcommand{\clo}{{\operatorname{closure}}}
\newcommand{\dist}{\operatorname{dist}}
\newcommand{\B}{\mathcal{B}}    % model
\newcommand{\calL}{\mathcal{L}} % LTI model class
\newcommand{\io}{\text{\rm i/o}}  
\newcommand{\tin}{\text{\rm i}}   % input
\newcommand{\tout}{\text{\rm o}}  % ouput
\newcommand{\ttm}{{\tt m}} % # inputs
\newcommand{\ttp}{{\tt p}} % # outputs
\newcommand{\ttl}{{\tt l}} % # latent variables
\newcommand{\ttd}{{\tt d}} % # latent variables
\newcommand{\ls}{\text{\rm ls}}
\newcommand{\wls}{\text{\rm wls}}
\newcommand{\als}{\text{\rm als}}
\newcommand{\tls}{\text{\rm tls}}
\newcommand{\gtls}{\text{\rm gtls}}
\newcommand{\wtls}{\text{\rm wtls}}
\newcommand{\fwtls}{\text{\rm fwtls}}
\newcommand{\rml}{\text{\rm l}}   % left, ch2
\newcommand{\rmr}{\text{\rm r}}   % right, ch2
\newcommand{\matlab}{{\sc Matlab}}
%--------------------------------------------------------------
\theoremstyle{remark}
\newtheorem{note}{Note} 
%--------------------------------------------------------------
\newcommand{\ie}{{\em i.e.}}
\newcommand{\eg}{{\em e.g.}}
\newcommand{\cf}{{\em cf.}}
\newcommand{\etc}{{\em etc}}
\newcommand{\etal}{{\em et~al.}}
%--------------------------------------------------------------
\renewcommand{\baselinestretch}{1}
\topmargin=-2.2cm  \oddsidemargin=-1cm \evensidemargin=-1cm 
\textheight=25.8cm   \textwidth=18cm
%--------------------------------------------------------------
\newcommand{\funinput}[1]{\medskip{{\small\tt\underline{#1}}}\verbatiminput{#1}}
\newcommand{\important}[1]{\begin{center}\framebox[.95\textwidth]{\begin{minipage}{.9\textwidth}\smallskip#1\smallskip\end{minipage}}\end{center}}
%--------------------------------------------------------------

\begin{document}

\title{A \matlab\ toolbox for weighted total least squares approximation}
\author{Ivan Markovsky and Sabine Van Huffel}
\maketitle
\date

\begin{abstract}
The toolbox solves a variety of approximate modeling problems for linear static models. The model can be parameterized in kernel, image, or input/output form and the approximation criterion, called misfit, is a weighted norm between the given data and data that is consistent with the model. There are three main classes of functions in the toolbox: transformation functions, misfit computation functions, and approximation functions. The approximation functions derive an approximate model from data, the misfit computation functions are used for validation and comparison of models, and the transformation functions are used for deriving one model representation from another.
\end{abstract}

\section{Introduction}\la{s1}

\subsection*{Purpose}

The weighted total least squares toolbox contains \matlab\ functions (m-files) for data approximation by linear static models. The data is a collection of~$N$, $\ttd$-dimensional real vectors $d_1,\ldots,d_N\in\R^\ttd$, gathered in a matrix $D:=\bmx d_1 & \cdots & d_N \emx\in\R^{\ttd\times N}$, and a linear static model~$\B$ for~$D$ is a subspace of~$\R^\ttd$. The natural number $\ttm:=\dim(\B)$ is a measure of the model complexity and $\calL^\ttd_{\ttm,0}$ denotes the set of all linear static models with $\ttd$ variables of dimension {\em at most\/}~$\ttm$.

A linear static model $\B\in\calL^\ttd_{\ttm,0}$ can be represented as a kernel or image of a matrix or in what is called an input/output form, see Table~\ref{t1}.
\begin{table}[htb!]
\centering
\caption{Representations of linear static models.} \la{t1}
\vskip.15cm
\begin{tabular}{|l|cl|l|}
\hline
Representation & \multicolumn{2}{|c|}{Parameter} & Model corresponding to that representation and parameter \\
\hline
Kernel         & $R\in\R^{g\times\ttd}$    &nonunique & $\B = \ker(R)   = \{\, d\in\R^\ttd \ | \ Rd = 0 \,\}$ \\
Image          & $P\in\R^{\ttd\times\ttl}$ &nonunique & $\B = \image(P) = \{\, d\in\R^\ttd \ | \ d = Pl, \ l\in\R^\ttl \,\}$ \\
Input/output   & $X\in\R^{\ttm\times\ttp}$ &unique     & $\B = \B_\io(X) = \{\, d = \col(d_\tin,d_\tout)\in\R^\ttd \ | \ d_\tout = X^\top d_\tin , \ d_\tin\in\R^\ttm \,\}$ \\
\hline
\end{tabular}
\end{table}

A representation of the model yields a parameterization. The model is described by equations that depend on parameter and to a given parameter corresponds a unique model. For a given model and a chosen representation, however, the corresponding parameter might not be unique. The nonuniqueness of the kernel and image representations are due to nonminimality of the representation ($R$ being not full row rank and $P$ being not full column rank) and due to freedom in the choice and of basis ($\ker(R) = \ker(UR)$ and $\image(P) = \image(PV)$ for any nonsingular matrices $U\in\R^{g\times g}$ and $V\in\R^{\ttl\times\ttl}$). The nonuniqueness in an input/output representation is due to the various possible partitionings of the variables into inputs (free variables) and outputs (dependent variables). The representation $\B_\io(X)$, however, has a fixed input/output partition, so that the parameter~$X$ is unique.

We use the short hand notation
$$
D = \bmx d_1 & \cdots & d_N\emx \in \B \quad :\iff \quad d_i\in\B, \quad\text{for } i=1,\ldots,N.
$$
If $D\in\B$, the model~$\B$ fits the data~$D$ exactly. If $D\not\in\B$, the model~$\B$ fit the data~$D$ only approximately. For optimal approximate modeling a distance function, called {\em misfit\/}, is defined as follows:
\begin{equation}
M_\wtls\big(\bmx d_1 & \cdots & d_N\emx,\B\big) := \min_{\hat d_1,\ldots,\hat d_N\in\B} \  \sqrt{ \sum_{i=1}^{N} (d_i - \hat d_i)^{\top} W_i (d_i - \hat d_i)} , \la{Mwtls}\tag{Mwtls}
\end{equation}
where $W_1,\ldots,W_N$ are given positive definite matrices. The weighted total least squares (WTLS) misfit $M_\wtls(D,\B)$ between the data~$D$ and a model $\B\in\calL_{\ttm,0}^\ttd$ is a measure of how much the model fails to fit the data exactly. The considered approximate modeling problem is:
\important{Given the data matrix $D = \bmx d_1 & \cdots & d_N\emx\in\R^{\ttd\times N}$, a complexity bound~$\ttm$, and weight matrices  $W_1,\ldots,W_N$, find an approximate model
\begin{equation}
\hat\B_{\wtls} := \arg\min_{\hat\B\in\calL_{\ttm,0}^\ttd} M_\wtls(D,\hat\B).  \la{wtls}\tag{WTLS}
\end{equation}\vskip-.25cm}
The special cases listed in Table~\ref{t2} allow for special solution methods and are treated separately.
\begin{table}[htb!]
\centering
\caption{Special cases of the weighted total least squares problem~\re{wtls}.} \la{t2}
\vskip.15cm
\begin{tabular}{|ll|l|l|}
\hline
Special case & & Name & Acronym \\
\hline
$W_i = \sigma^2I$        & $\sigma\in\R_+$                & total least squares & TLS \\
$W_i = \diag(w)$         & $w\in\R_+^\ttd$                & element-wise generalized total least squares & EWGTLS \\
$W_i = W$                & $W>0$                          & generalized total least squares & GTLS \\
$M_{\wtls} = M_{\gtls2}$ & $W_\rml,W_\rmr>0$ diagonal     & EWGTLS with two side weighting & EWGTLS2 \\
$M_{\wtls} = M_{\gtls2}$ & $W_\rml,W_\rmr>0$              & GTLS with two side weighting & GTLS2 \\
$W_i = \diag(w_i)$       & $w_i\in\R_+^\ttd$              & element-wise weighted total least squares & EWTLS \\
\hline
\end{tabular}
\begin{equation*}
M_{\gtls2}(D,\B) := \min_{\hat D\in\B}\norm{\sqrt{W_\rml}(D-\hat D)\sqrt{W_\rmr}}_\F %\la{gtls2}\tag{GTLS2}
\end{equation*}
\end{table}

\begin{note}
The following weighted total least squares problem, called fully weighted total least square (FWTLS) problem
$$
\hat\B_{\fwtls} := \arg\min_{\hat\B\in\calL_{\ttm,0}^\ttd} \min_{\hat D\in\hat\B} \mtov^\top(D-\hat D) W \mtov(D-\hat D),\quad \text{where}\quad W\in\R^{\ttd N\times\ttd N}
$$
is also considered. It includes~\re{wtls} as a special case with $W=\blkdiag(W_1,\ldots,W_N)$. The FWTLS problem, however, does not allow for efficient computational methods and its solution is prohibitive already for small sample size problems (say $\ttd=10$ and $N=100$). For this reason the FWTLS problem is not the central problem of interest and is included only for completeness. 
\end{note}

\subsection*{Algorithms}

The special cases listed in Table~\ref{t2} have increased generality from top to bottom. The more general the problem is, however, the more computationally expensive its solution is. The TLS and GTLS problems allow for analytic solutions in terms of the singular value decomposition (SVD). The more general EWTLS, WTLS, and FWTLS problems (currently) have no such analytic solutions and use less robust iterative solution methods. 

The SVD method is computationally faster than the alternative iterative optimization methods and theoretically characterizes all globally optimal solutions. In contrast, the iterative optimization methods (used in the package) compute one locally optimal solution. (The algorithm of Premoli--Rastello \cite{PR01,MRPKV02} is not globally convergent to a local solution, so that for particular initial approximations this method might not converge to a local solution. In such cases the algorithm diverges or oscillates.)

The GTLS-type problems (EWGTLS, GTLS, EWGTLS2, and GTLS2) are solved in the package via a transformation technique. The data matrix is appropriately scaled and the corresponding TLS problem is solved for the scaled data. Then the solution of the original problem is recovered from the solution of the transformed TLS problem via the inverse transformation. 

The general WTLS problem is solved via local optimization methods. The following algorithms are used/implemented in the package:
\begin{enumerate}
\item classical local optimization methods (from the Optimization Toolbox of \matlab),
\item an alternating least squares algorithm,
\item the algorithm of Premoli--Rastello.
\end{enumerate}

\subsection*{Implementation}

The implementation is in \matlab\ code. For the problems with analytic solution, the \matlab\ code is expected to compute a solution nearly as fast as an alternative code in C or FORTRAN. The general WTLS algorithms, however, are expected to benefit in terms of execution time if implemented in C or FORTRAN.