Mercurial > hg > octave-lojdl
changeset 7681:b1c1133641ee
Add the spaugment function
author | David Bateman <dbateman@free.fr> |
---|---|
date | Wed, 02 Apr 2008 14:08:28 -0400 |
parents | a0ec02774303 |
children | f7474c83c782 |
files | doc/ChangeLog doc/interpreter/sparse.txi scripts/ChangeLog scripts/sparse/Makefile.in scripts/sparse/spaugment.m |
diffstat | 5 files changed, 112 insertions(+), 3 deletions(-) [+] |
line wrap: on
line diff
--- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -1,3 +1,7 @@ +2008-04-02 David Bateman <dbateman@free.fr> + + * interpreter/sparse.txi: Document spaugment. + 2008-03-26 Rafael Laboissiere <rafael@debian.org> * interpreter/mkoctfile.1: Remove spurious whitespace before macros
--- a/doc/interpreter/sparse.txi +++ b/doc/interpreter/sparse.txi @@ -474,7 +474,7 @@ @item Linear algebra: @dfn{matrix_type}, @dfn{normest}, @dfn{condest}, @dfn{sprank} -@c @dfn{spaugment} + @dfn{spaugment} @c @dfn{eigs}, @dfn{svds} but these are in octave-forge for now @item Iterative techniques: @@ -825,6 +825,11 @@ @DOCSTRING(symbfact) +For non square matrices, the user can also utilize the @code{spaugment} +function to find a least squares solution to a linear equation. + +@DOCSTRING(spaugment) + @node Iterative Techniques, Real Life Example, Sparse Linear Algebra, Sparse Matrices @section Iterative Techniques applied to sparse matrices
--- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -1,5 +1,8 @@ 2008-04-02 David Bateman <dbateman@free.fr> + * sparse/spaugment.m: New function + * sparse/Makefile.in (SOURCES): Add it here. + * plot/__gnuplot_ginput__.m: Use the gnuplot stream itself for communication rather than a chat file if mkfifo is not available. * plot/gnuplot_drawnow.m: Open stream with popen2 to allow two way
--- a/scripts/sparse/Makefile.in +++ b/scripts/sparse/Makefile.in @@ -33,8 +33,8 @@ INSTALL_DATA = @INSTALL_DATA@ SOURCES = colperm.m etreeplot.m gplot.m nonzeros.m normest.m \ - pcg.m pcr.m spalloc.m spconvert.m spdiags.m speye.m spfun.m \ - sphcat.m spones.m sprand.m sprandn.m sprandsym.m spstats.m \ + pcg.m pcr.m spalloc.m spaugment.m spconvert.m spdiags.m speye.m \ + spfun.m sphcat.m spones.m sprand.m sprandn.m sprandsym.m spstats.m \ spvcat.m spy.m treeplot.m DISTFILES = $(addprefix $(srcdir)/, Makefile.in $(SOURCES))
new file mode 100644 --- /dev/null +++ b/scripts/sparse/spaugment.m @@ -0,0 +1,97 @@ +## Copyright (C) 2008 David Bateman +## +## This file is part of Octave. +## +## Octave is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or (at +## your option) any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{s} =} spaugment (@var{a}, @var{c}) +## Creates the augmented matrix of @var{a}. This is given by +## +## @example +## [@var{c} * eye(@var{m}, @var{m}),@var{a}; @var{a}', zeros(@var{n}, +## @var{n})] +## @end example +## +## @noindent +## This is related to the leasted squared solution of +## @code{@var{a} \\ @var{b}}, by +## +## @example +## @var{s} * [ @var{r} / @var{c}; x] = [@var{b}, zeros(@var{n}, +## columns(@var{b})] +## @end example +## +## @noindent +## where @var{r} is the residual error +## +## @example +## @var{r} = @var{b} - @var{a} * @var{x} +## @end example +## +## As the matrix @var{s} is symmetric indefinite it can be factorized +## with @code{lu}, and the minimum norm solution can therefore be found +## without the need for a @code{qr} factorization. As the residual +## error will be @code{zeros (@var{m}, @var{m})} for under determined +## problems, and example can be +## +## @example +## @group +## m = 11; n = 10; mn = max(m ,n); +## a = spdiags ([ones(mn,1), 10*ones(mn,1), -ones(mn,1)],[-1,0,1], m, n); +## x0 = a \ ones (m,1); +## s = spaugment (a); +## [L, U, P, Q] = lu (s); +## x1 = Q * (U \ (L \ (P * [ones(m,1); zeros(n,1)]))); +## x1 = x1(end - n + 1 : end); +## @end group +## @end example +## +## To find the solution of an overdetermined problem needs an estimate +## of the residual error @var{r} and so it is more complex to formulate +## a minimum norm solution using the @code{spaugment} function. +## +## In general the left division operator is more stable and faster than +## using the @code{spaugment} function. +## @end deftypefn + +function s = spaugment (a, c) + if (nargin < 2) + if (issparse (a)) + c = max (max (abs (a))) / 1000; + else + if (ndims (a) != 2) + error ("spaugment: expecting 2-dimenisional matrix") + else + c = max (abs (a(:))) / 1000; + endif + endif + elseif (!isscalar (c)) + error ("spaugment: c must be a scalar"); + endif + + [m, n] = size (a); + s = [ c * speye(m, m), a; a', sparse(n, n)]; +endfunction + +%!test +%! m = 11; n = 10; mn = max(m ,n); +%! a = spdiags ([ones(mn,1), 10*ones(mn,1), -ones(mn,1)],[-1,0,1], m, n); +%! x0 = a \ ones (m,1); +%! s = spaugment (a); +%! [L, U, P, Q] = lu (s); +%! x1 = Q * (U \ (L \ (P * [ones(m,1); zeros(n,1)]))); +%! x1 = x1(end - n + 1 : end); +%! assert (x1, x0, 1e-10)