INTERNATIONAL WORKSHOP ON
NUMERICAL LINEAR ALGEBRA
RECENT TRENDS IN
RECENT TRENDS IN
INVITED TALKS
Abstract.
We report the results of some recent joint researches with L. Gemignani,
V. Pan and F. Tisseur where semiseparable matrices play an important role.
A first result concerns a method for the solution of the generalized eigenvalue
problem (A-xB)v=0, where A is a real symmetric tridiagonal matrix and B
is a real (nondefinite) diagonal matrix, based on Aberth's iteration where
the evaluation of the ratio p(x)/p'(x), p(x)=det(A-xB), is performed
in a stable and robust way by exploiting the QR factorization of A-xB and
the semiseparability of Q. A second result concerns the QR iteration applied
to an nxn companion matrix F or to a matrix of the kind D+uv^T, where D
is an nxn real diagonal matrix and $u,v\in R^n$. We show that the QR iteration,
either with linear or with quadratic shift, generates a sequence of semiseparablematrices
so that each QR step can be implemented with O(n) ops. Applications to
approximating zeros of polynomials are shown.
Abstract.
Tikhonov regularization is one of the most popular regularization techniques
for the solution of linear ill-posed problems.
Large-scale ill-posed problems have to be solved by iterative techniques,
and we review available and present new iterative
methods. We consider both the cases when the norm of the error in the
available observations is known and when it is not. A rational for using
nonsingular regularization operators is also discussed.
Abstract.
Block p-cyclic matrices arise in a wide variety of applications. This
includes the solution of the approximation to elliptic
partial differential equations and computing the stationary probability
distribution of Markov matrices. In this talk, we
shall show how we can specialize Krylov space methods to this class
of matrices and show the simplifications that can occur. A
substantial savings in iterations can be made using these devices.The
methods proposed are particularly useful in the context of
Domain Decomposition. Numerical experiences will be reported.
Abstract.
The unitary Hessenberg QR (uhqr) algorithm is fundamental for statistical
signal analysis, as is the closely related inverse
algorithm (iuhqr). These algorithms are analogous with algorithms,
tqt and itqr, for real symmetric tridiagonal matrices.
No claim of numerical stability for uhqr was made when it was introduced
in 1986. Indeed, an open problem was to make uhqr perform as well
as tqr. Here we introduce a device which appears, on the basis of
a large number of numerical experiments, to solve the problem. Indeed,
Michael Stewart has proved this in detail for a slight variation of our
algorithm.
Abstract.
We consider the problem of calculating the eigenvalues and singular
values of complex Hankel matrices, and we present fast
O(n^2log(n)) algorithms for the computation. Our techniques are based
on a generalized Lanczos iteration and the Fast Fourier
Transform.
Abstract.
It is well-known that the proper way to deal with perturbations of eigenvectors
and singular vectors is to bound the distance between invariant or singular
subspaces. From this point of view having close subspaces is equivalent
to the existence of close bases. We show that this is no longer true for
simultaneous bases of singular subspaces, i.e. for pairs of orthonormal
bases of
corresponding left and right singular subspaces appearing in a same
singular value decomposition. A new perturbation theory for simultaneous
singular bases is developed, both in the absolute and the relative setting,
showing that important differences with respect to the sensitivity of singular
subpaces appear only in the absolute case. These results guarantee, in
particular, that,
unlike classical algorithms as QR or divide and conquer, high relative
accuracy algorithms for the singular value decomposition
producing multiplicative backward errors do compute reliably such simultaneous
bases of singular subspaces. Illustrative numerical experiments are presented.
Abstract.
In this paper, we focus on image deconvolution and image reconstruction
problems where a sought image is recovered from
degraded observed data. The solution is defined to be the minimizer
of an objective function combining a data-fidelity term
and a edge-preserving, convex regularization term. To accelerate the
computation of the estimate, two forms of half-quadratic
regularization, multiplicative and additive, are often used. The goal
of this paper is to compare both theoretically and
experimentally the effectiveness and the efficiency of these two forms.
We provide a theoretical and experimental analysis of the
speed of convergence that they allow. We also propose a method applying
pertinent preconditioning to an adapted half-quadratic equivalent form
of the objective function. We focus specifically on Huber regularization.
We exhibit the possibility get very fast calculations while preserving
the edges in the solution. Finally, we study the method for the MRI image
reconstruction and other constrained regularizations in computer vision
problems.
Abstract.
In this talk we provide a unified approach to derive fast algorithms
to compute discrete Cosine and Sine transforms I - IV.
For the Sine I/II and Cosine I/II transforms a unified derivation
based on the concept of the comrade matrix is presented. The
comrade matrices associated with different versions of the transforms
differ in only a few boundary elements; hence, in each
case algorithms can be derived in a unified manner. The algorithm is
then modified to compute Sine III/IV and Cosine III/IV
transforms as well. The resulting algorithms for the versions III/IV
are direct and recursive, such algorithms were missing in
the existing literature. Finally, formulas reducing Cosine and
Sine transforms of the types III and IV to each other are
presented.
Abstract.
Many discretized large-scale ill-posed problems yield matrices that
are nonsymmetric and possibly nonsquare. Traditionally, the
solution methods of choice for such problems have been iterative methods
applied to the normal equations or a modification of the normal equations
obtained by Tikhonov regularization; the problems solved numerically are
large linear systems of equations with a symmetric matrix. We have found
that by respecting the nonsymmetry of the problem, i.e., by solving the
normal equations and instead applying iterative methods, such as the GMRES,
BiCG or QMR methods, to the nonsymmetric linear system one may obtain faster
convergence and an approximate solution of higher quality than when solving
the normal equations by the CG method. In fact, in applications to image
restoration the new methods can be competitive with TV-norm based methods.
Abstract.
In this talk, we present an inverse free Krylov subspace method for
finding some extreme eigenvalues of the symmetric definite
generalized eigenvalue problem Ax = lambda B x. The basic
method takes a form of inner-outer iterations and involves no
inversion of B or any shift-and-invert matrix. A convergence analysis
is presented that leads to a preconditioning scheme for
accelerating convergence through some equivalent transformations of
the eigenvalue problem. Numerical examples are given to
illustrate the convergence properties and to demonstrate the competitiveness
of the method.
SHORT COMMUNICATIONS
Abstract.
In this work we present a parallel algorithm for solving a system of
linear equations
Tx=b, (1)
where $T\in \mathbb{R}^{n\times n}$ is a real symmetric Toeplitz matrix.
There exist fast algorithms for solving Toeplitz linear systems which are based on the \emph{displacement rank} property of these matrices. But the dependency among very low grain operations makes them difficult to parallelize efficiently. If we use distributed memory architectures, those parallel algorithms have large communication costs with respect to their small computation costs. The communications are mainly due to the form of the displacement matrices involved in the displacement representation of a Toeplitz matrix and its influence in the algorithm.
In order to overcome the high cost of the communications mentioned above,
we propose to transform the Toeplitz matrix into a
generalized Cauchy matrix. The generalized Cauchy matrices are also
structured matrices. The displacement matrices involved in their displacement
representation are diagonal. This fact let us avoid a lot of communications
in the parallel algorithm. Furthermore, the Cauchy-like matrices arising
from the transformation of a symmetric Toeplitz matrix have a considerable
sparsity that can be exploited.
On the other hand, our main goal has been to develop a parallel algorithm
to solve (1) efficiently in general purpose distributed
parallel environments. Also, we are interested in using standard and
public domain libraries in order to get suitable, efficient
and portable codes. We have used LAPACK and ScaLAPACK libraries based
on the message-passing library MPI. Besides, the availability of MPI versions
for shared memory architectures allows to run our codes on multiprocessor
boards. We will show how to use the routines and the distribution of data
based on the ScaLAPACK logical topology in order to get an efficient and
portable algorithm.
We will also show in the exposition how to apply the parallel algorithm
to solve an inverse filtering of multichannel systems
problem arising in digital processing of sound signals. This interesting
problem is mathematically modeled by the minimization
of a squared error expression involving Toeplitz matrices. It can be
approached by solving a system of linear equations where the system matrix
is Toeplitz-block with symmetric Toeplitz blocks. With a few modifications
of the basic parallel algorithm it is
possible to use it in the Toeplitz-block case.
Abstract.
We present analyticity results and new algorithms for computing power
series approximations to the function A^{-1}(t)b(t) for
$b: R \rightarrow R^n$ and $A: R \rightarrow R^{n\times n}$. These
procedures have application whenever systems A(t)x=b(t) need be solved
for $t=t_0...t_p$ and p is large. We demonstrate these new methods to be
competitive with both direct factorizations and preconditioned iterative
solvers. The bulk of our computation can be done offline and hence
these methods are particularly suitable for real time applications. In
the special case that A is a polynomial, $A(t)=A_0+...+A_pt^p$, there is
a particularly fast algorithm. These matrix polynomials arise frequently
such as in solving discretized linear ode's in Fourier space, for example
frequency response equations. We describe current research directions into
the more general case where $A: R^m \rightarrow R^{n\times n}$.
Abstract.
In this talk we discuss a generalization of singular value decomposition
to the 3 and more dimensions. If matrices of
singular vectors have full column rank then it is possible to construct
an algorithm for this decomposition with O(sr+r^5) arithmetical operations,
where r is the rank of this decomposition and s the product of all dimensions
in the problem.
In this talk we will discuss the stability of this approach and present
some new ideas how to reduce the arithmetical complexity
of this method.
Abstract.
The computation of the Singular Value Decomposition (SVD) of structured
matrices such as Toeplitz, Hankel, Cauchy or
Vandermonde matrices has become an important line of research in numerical
linear algebra. In this work the problem of inversion in the context of
the computation of curve intersections is considered. Although this problem
has usually been dealt with in the field of exact rational computations
and in that case it can be solved by using Gaussian elimination, when one
has to work in finite precision arithmetic the problem leads to the computation
of the SVD of a Sylvester matrix, a different type of structured matrix
widely used in computer algebra. In addtition only a small part of the
SVD is needed, which shows the interest of having special algorithms for
this situation.
Abstract.
In [1] a method for localizating the real parts of the eigenvalues of
a real matrix was proposed. This method can be considered as an alternative
to the Gerschgorin circles and in [2] it has been proved that it also provides
an alternative to the ovals of
Cassini. In this communication we compare the accuracy of the Gerschgorin
circles and our alternative method for some special
classes of matrices and we present some applications obtained by combining
both localization methods.
[1] Peña, J.M: A class of P-matrices with applications to the
localization of the eigenvalues of a real matrix, SIAM Journal on
Matrix Analysis and its Applications 22(2001), 1027-1037.
[2] Peña, J.M.: On an alternative to Gerschgorin circles and
ovals of Cassini. To appear in Numerische Mathematik.
Abstract.
A symplectic Arnoldi method for a large and sparse matrix eigenvalue
problem is presented. The method is based on a
significant updating of both Gram-Schmidt and Arnoldi process. The
motivation lies in constructing structure-preserving Krylov
subspace methods for a large and sparse Hamiltonian matrix. Such methods
are appropriate for the computation of several Hamiltonian matrix eigenvalues
and eigenvectors or when low-rank approximations to the solutions of a
continuous-time algebraic Riccati equations are required. A general framework,
based on subspace Krylov methods is given. New algorithms such symplectic
Gram-Schmidt, modified symplectic Gram-Schmidt, symplectic Arnoldi, modified
symplectic Arnoldi and Hamiltonian symplectic Arnoldi are obtained. Some
of their properties are established. Some of their potentialities and differences
with some existing methods are mentioned. The symplectic Lanczos algorithm
is obtained from Hamiltonian symplectic Arnoldi as a particular case.
Abstract.
We construct dynamical systems which evolution converge to the eigenvectors
(eigenvalues) of a general square matrix, not
neccesarily symmetric. The properties of the systems are analyzed and
some applications are presented.