function NearestOrthFun()
Introduction
Suppose $A$ is a given $m \times n$ matrix. The problem of finding the orthonormal matrix $Q$ nearest to $A$ is well-known [4] (by an orthonormal matrix, we mean a matrix with orthonormal columns). Specifically, it is the following problem: $$ Q = \mbox{argmin}_W ||A-W||_\mbox{fro} \quad\mbox{subject to}\quad W^T W = I. $$ This problem is a special case of the orthogonal Procrustes problem [3, p. 327]. There are different ways to find the unique solution to this problem. One way that expresses $Q$ explicitly requires finding a matrix inverse square root [6] as follows:
$$ Q = A(A^T A)^{-1/2}. $$
However, the simplest (and possibly a numerically more stable) way is to compute the singular value decomposition of $A$ and to replace the singular values by ones. In other words, if $A = U S V^T$, then $Q = U V^T$ is the nearest orthonormal matrix to $A$. This is actually the unitary factor in the polar decomposition of $A$ [4].
Our goal is to generalize this situation to the continuous case. Suppose that a set of functions defined on a domain are given. We know how to use Chebfun's qr
command to compute orthogonal functions (stored as the columns of a qusimatrix Q2
) corresponding to our given functions. See e.g. [2, Ch. 6] for details. Here, we are looking for the nearest set of orthonormal functions. In Chebfun this is also straightforward. We just need to create a quasimatrix out of our given set of functions and use Chebfun's overloaded svd
command [1]. The function [Q,Q2] = nearestOrtho(A)
in the following computes both Q
and Q2
.
The first six monomials
A = cheb.gallery('vandermonde'); % A = [1, x, x.^2, x.^3, x.^4, x.^5]; [Q,Q2] = nearestOrtho(A);
Departure from orthogonality in the columns of A = 1.35 Departure from orthogonality in the columns of Q = 2.1e-15 Departure from orthogonality in the columns of Q2 = 1.7e-15 The distance between A and Q2 = 1.92 The distance between A and its closest orthonormal quasimatrix = 1.69
Here, we see Legendre polynomials scaled to have norm 1 (left) and the nearest set of orthonormal functions (right). The first polynomial among these nearest orthonormal polynomials is approximately the following: $ p(x) = 0.9200 - 0.9035 x^2 + 0.3084 x^4$.
Chebyshev-Vandermonde quasimatrix on the interval [0, 1]
A = cheb.gallery('vandercheb'); A = restrict(A,[0,1]); [Q,Q2] = nearestOrtho(A);
Departure from orthogonality in the columns of A = 1.00 Departure from orthogonality in the columns of Q = 2.8e-15 Departure from orthogonality in the columns of Q2 = 1.7e-15 The distance between A and Q2 = 2.48 The distance between A and its closest orthonormal quasimatrix = 1.62
Mixed algebraic and trigonometric example
x = chebfun('x'); A = [1, cos(x), sin(x.^2), x.^3, x.^4, x.^5]; [Q,Q2] = nearestOrtho(A);
Departure from orthogonality in the columns of A = 2.69 Departure from orthogonality in the columns of Q = 5.1e-15 Departure from orthogonality in the columns of Q2 = 5.7e-15 The distance between A and Q2 = 2.43 The distance between A and its closest orthonormal quasimatrix = 1.95
Non-smooth functions on the interval [0,10]
A = [cheb.gallery('stegosaurus'), cheb.gallery('wiggly'), cheb.gallery('blasius')]; [Q,Q2] = nearestOrtho(A);
Departure from orthogonality in the columns of A = 194.38 Departure from orthogonality in the columns of Q = 1.4e-15 Departure from orthogonality in the columns of Q2 = 3.8e-16 The distance between A and Q2 = 13.63 The distance between A and its closest orthonormal quasimatrix = 13.23
A further comment
We mentioned in the introduction that the nearest orthonormal matrix $Q$ to a given matrix $A$ is the orthonormal factor of a polar deomposition $A = QH$. The other factor $H$ is a Hermitian matrix. Analogously, in the case considered in this Example of functions rather than vectors, we are computing the nearest orthognormal quasimatrix $Q$ to a given quasimatrix $A$, and that could be regarded again as half of a polar decomposition $A = QH$. Here both $A$ and $Q$ are quasimatrices, but $H$ is again a fully discrete matrix. The Chebfun team may introduce at a later date a Chebfun command POLDEC to compute the polar decomposition (analogous to Higham's POLDEC command for matrices in the Matrix Computations Toolbox [5]), in which case, this example could be simplified.
function [Q,Q2] = nearestOrtho(A) [U,S,V] = svd(A); Q = U*V'; % nearest quasimatrix to A with orthonormal columns [Q2,R] = qr(A); m = size(A,2); figure; subplot(1,2,1); plot(Q2); v=axis; v(4) = -v(3); axis(v); % make the scales look better grid on; title('QR orthonormalization') subplot(1,2,2); plot(Q); axis(v); grid on; title('Optimal orthonormalization') fprintf(['Departure from orthogonality in the columns' ... ' of A = %3.2f\n'], norm(A'*A-eye(m))) fprintf(['Departure from orthogonality in the columns' ... ' of Q = %3.2g\n'], norm(Q'*Q-eye(m))) fprintf(['Departure from orthogonality in the columns' ... ' of Q2 = %3.2g\n'], norm(Q2'*Q2-eye(m))) fprintf('The distance between A and Q2 = %3.2f\n', norm(A-Q2,'fro')) fprintf(['The distance between A and its closest orthonormal' ... ' quasimatrix = %3.2f\n'], norm(A-Q,'fro')) end
References
-
Z. Battles and L.N. Trefethen, An extension of MATLAB to continuous functions and operators, SIAM J. Scientific Computing 25 (2004) 1743-1770.
-
T.A. Driscoll, N. Hale, and L.N. Trefethen, editors, Chebfun Guide, Pafnuty Publications, Oxford, 2014.
-
G.H. Golub and C.F. Van Loan, Matrix Computations, 4th edition, The Johns Hopkins University Press, 2013.
-
N.J. Higham, Matrix nearness problems and applications, pp. 1-27, in Applications of Matrix Theory, Inst. Math. Appl. Conf. Ser. New Ser., Vol. 22, Oxford Univ. Press, New York, 1989.
-
N.J. Higham, The matrix computation toolbox, http://www.ma.man.ac.uk/~higham/mctoolbox/, 2008.
-
B.K.P. Horn, H.M. Hilden, and S. Negahdaripour, Closed-form solution of absolute orientation using orthonormal matrices, Journal of the Optical Society A, 5 (1988) 1127-1135.
end