juno-machine/src/RedundantLSolve.i3

 Copyright (C) 1992, Digital Equipment Corporation                         
 All rights reserved.                                                      
 See the file COPYRIGHT for a full description.                            
                                                                           
 Last modified on Thu May  2 09:36:00 PDT 1996 by heydon                   
      modified on Fri Dec  2 15:35:09 PST 1994 by gnelson                  
      modified on Mon Oct 31 18:13:52 PST 1994 by isard                    


INTERFACE RedundantLSolve;

 This interface provides a procedure P that approximately solves
   a set of simultaneous linear equations.  P is designed to be
   used as a subroutine to solve simultaneous nonlinear equations
   by Newton interation: at each stage of the iteration, the nonlinear
   equations are approximated by linear ones around the current point,
   these linear equations are solved, and the current point is moved
   to the solution found.


   The reason that P finds only an approximate solution instead
   of an exact one is that the approximate solution works better
   in the case that the nonlinear equations are redundant but consistent
   (for example, x^2=4 AND x^3=8).  In this case the linear approximation
   often produces inconsistent or very ill-conditioned linear systems.


   So, instead of solving



      x : Ax = b


    the procedure P solves



       x : | b - Ax | < |b| / 2


    That is, it sets x so that the norm of b - Ax is less than half the
    norm of b.  This criterion guarantees that progress in made in the
    overall Newton iteration on the non-linear problem.


    The procedure P attempts to find a small x that satisfies the
    constraint.  It is not guaranteed to find the smallest x in any
    precise sense, but if A is very ill-conditioned (as in the case
    where the nonlinear problem is redundant but consistent) the
    procedure avoids the directions of motion corresponding to the
    small singular values of A, and if A is very underconstrained,
    so that there are many dimensions worth of xs that will solve the
    constraint, then the Gram-Schmidt algorithm is used to find
    the one that is shortest in the Euclidean norm. 


FROM JunoValue IMPORT Real;

IMPORT Wr;

TYPE
  Vector = ARRAY OF Real;
  Matrix = ARRAY OF Vector;

PROCEDURE P(
  m, n: CARDINAL;
  VAR (*INOUT*) a: Matrix;
  VAR (*OUT*) x: Vector);

 Set x to the approximate solution of Ax=b, where matrix A is stored
   as a[0..m-1,0..n-1] and vector b is stored as a[0..m-1,n]. 

VAR logWr: Wr.T := NIL;

 If the internal variable RedundantLSolve.debug exceeds zero,
   information is written to logWr each time P is called. 


PROCEDURE SetGramSchmidt(on: BOOLEAN);

 By default, calls to P use the Gram-Schmidt algorithm to compute a
   minimal solution in the underconstrained case. This procedure can be used
   to disable (or re-enable) that process. 

END RedundantLSolve.





interface JunoValue is in: