Package gov.nih.mipav.model.algorithms

Class LSQR

  • All Implemented Interfaces:
    java.awt.event.ActionListener, java.awt.event.WindowListener, java.lang.Runnable, java.util.EventListener
    public class LSQR
extends AlgorithmBase

    • Constructor Detail

      • LSQR

        public LSQR()

    • Method Detail

      • runAlgorithm

        public void runAlgorithm()
        Description copied from class: AlgorithmBase
        Actually runs the algorithm. Implemented by inheriting algorithms.
        Specified by:
        runAlgorithm in class AlgorithmBase

      • aprod_ez

        public void aprod_ez​(LSQR.lsqr_solver_ez me,
                     int mode,
                     int m,
                     int n,
                     double[] x,
                     double[] y)

      • solve_ez

        public void solve_ez​(LSQR.lsqr_solver_ez me,
                     double[] b,
                     double damp,
                     double[] x,
                     int[] istop,
                     double[] se,
                     int[] itn,
                     double[] anorm,
                     double[] acond,
                     double[] rnorm,
                     double[] arnorm,
                     double[] xnorm)

      • LSQR

        public void LSQR​(LSQR.lsqr_solver_ez me,
                 int m,
                 int n,
                 double damp,
                 boolean wantse,
                 double[] u,
                 double[] v,
                 double[] w,
                 double[] x,
                 double[] se,
                 double atol,
                 double btol,
                 double conlim,
                 int itnlim,
                 java.io.RandomAccessFile raFile,
                 int[] istop,
                 int[] itn,
                 double[] anorm,
                 double[] acond,
                 double[] rnorm,
                 double[] arnorm,
                 double[] xnorm)
        !> ! LSQR finds a solution \(x\) to the following problems: ! ! 1. Unsymmetric equations -- solve: $$ \mathbf{A} \cdot \mathbf{x} = \mathbf{b} $$ ! 2. Linear least squares -- solve (in the least-squares sense): $$ \mathbf{A} \cdot \mathbf{x} = \mathbf{b} $$ ! 3. Damped least squares -- solve (in the least-squares sense): $$ \left[ \begin{array}{c} \mathbf{A}\\ damp \cdot \mathbf{I} \end{array} \right] \cdot \mathbf{x} = \left[ \begin{array}{c} \mathbf{b}\\ \mathbf{0} \end{array} \right] $$ ! ! where A is a matrix with m rows and n columns, b is an ! m-vector, and damp is a scalar. (All quantities are real.) ! The matrix A is intended to be large and sparse. It is accessed ! by means of subroutine calls to aprod. ! ! The rhs vector b is input via u, and subsequently overwritten. ! ! Note: LSQR uses an iterative method to approximate the solution. ! The number of iterations required to reach a certain accuracy ! depends strongly on the scaling of the problem. Poor scaling of ! the rows or columns of A should therefore be avoided where ! possible. ! ! For example, in problem 1 the solution is unaltered by ! row-scaling. If a row of A is very small or large compared to ! the other rows of A, the corresponding row of ( A b ) should be ! scaled up or down. ! ! In problems 1 and 2, the solution x is easily recovered ! following column-scaling. Unless better information is known, ! the nonzero columns of A should be scaled so that they all have ! the same Euclidean norm (e.g., 1.0). ! ! In problem 3, there is no freedom to re-scale if damp is ! nonzero. However, the value of damp should be assigned only ! after attention has been paid to the scaling of A. ! ! The parameter damp is intended to help regularize ! ill-conditioned systems, by preventing the true solution from ! being very large. Another aid to regularization is provided by ! the parameter acond, which may be used to terminate iterations ! before the computed solution becomes very large. ! ! Note that x is not an input parameter. ! If some initial estimate x0 is known and if damp = 0, ! one could proceed as follows: ! ! 1. Compute a residual vector \( r_0 = b - A \cdot x_0 \). ! 2. Use LSQR to solve the system \( A \cdot \Delta x = r_0 \). ! 3. Add the correction dx to obtain a final solution \( x = x_0 + \Delta x \). ! ! This requires that x0 be available before and after the call ! to LSQR. To judge the benefits, suppose LSQR takes k1 iterations ! to solve A*x = b and k2 iterations to solve A*dx = r0. ! If x0 is "good", norm(r0) will be smaller than norm(b). ! If the same stopping tolerances atol and btol are used for each ! system, k1 and k2 will be similar, but the final solution x0 + dx ! should be more accurate. The only way to reduce the total work ! is to use a larger stopping tolerance for the second system. ! If some value btol is suitable for A*x = b, the larger value ! btol*norm(b)/norm(r0) should be suitable for A*dx = r0. ! ! Preconditioning is another way to reduce the number of iterations. ! If it is possible to solve a related system M*x = b efficiently, ! where M approximates A in some helpful way ! (e.g. M - A has low rank or its elements are small relative to ! those of A), LSQR may converge more rapidly on the system ! A*M(inverse)*z = b, ! after which x can be recovered by solving M*x = z. ! ! NOTE: If A is symmetric, LSQR should not be used! ! Alternatives are the symmetric conjugate-gradient method (cg) ! and/or SYMMLQ. ! SYMMLQ is an implementation of symmetric cg that applies to ! any symmetric A and will converge more rapidly than LSQR. ! If A is positive definite, there are other implementations of ! symmetric cg that require slightly less work per iteration ! than SYMMLQ (but will take the same number of iterations). ! !### Notation ! ! The following quantities are used in discussing the subroutine ! parameters: ! ! ! Abar = ( A ), bbar = ( b ) ! ( damp*I ) ( 0 ) ! ! r = b - A*x, rbar = bbar - Abar*x ! ! rnorm = sqrt( norm(r)**2 + damp**2 * norm(x)**2 ) ! = norm( rbar ) ! ! relpr = the relative precision of floating-point arithmetic ! on the machine being used. On most machines, ! relpr is about 1.0e-7 and 1.0d-16 in single and double ! precision respectively. ! ! ! LSQR minimizes the function rnorm with respect to x. ! !### References ! ! * C.C. Paige and M.A. Saunders, LSQR: An algorithm for sparse ! linear equations and sparse least squares, ! ACM Transactions on Mathematical Software 8, 1 (March 1982), ! pp. 43-71. ! ! * C.C. Paige and M.A. Saunders, Algorithm 583, LSQR: Sparse ! linear equations and least-squares problems, ! ACM Transactions on Mathematical Software 8, 2 (June 1982), ! pp. 195-209. ! ! * C.L. Lawson, R.J. Hanson, D.R. Kincaid and F.T. Krogh, ! Basic linear algebra subprograms for Fortran usage, ! ACM Transactions on Mathematical Software 5, 3 (Sept 1979), ! pp. 308-323 and 324-325. ! !### LSQR development ! ! * 22 Feb 1982: LSQR sent to ACM TOMS to become Algorithm 583. ! * 15 Sep 1985: Final F66 version. LSQR sent to "misc" in netlib. ! * 13 Oct 1987: Bug (Robert Davies, DSIR). Have to delete ! if ( (one + dabs(t)) <= one ) GO TO 200 ! from loop 200. The test was an attempt to reduce ! underflows, but caused w(i) not to be updated. ! * 17 Mar 1989: First F77 version. ! * 04 May 1989: Bug (David Gay, AT&T). When the second beta is zero, ! rnorm = 0 and ! test2 = arnorm / (anorm * rnorm) overflows. ! Fixed by testing for rnorm = 0. ! * 05 May 1989: Sent to "misc" in netlib. ! * 14 Mar 1990: Bug (John Tomlin via IBM OSL testing). ! Setting rhbar2 = rhobar**2 + dampsq can give zero ! if rhobar underflows and damp = 0. ! Fixed by testing for damp = 0 specially. ! * 15 Mar 1990: Converted to lower case. ! * 21 Mar 1990: d2norm introduced to avoid overflow in numerous ! items like c = sqrt( a**2 + b**2 ). ! * 04 Sep 1991: wantse added as an argument to LSQR, to make ! standard errors optional. This saves storage and ! time when se(*) is not wanted. ! * 13 Feb 1992: istop now returns a value in [1,5], not [1,7]. ! 1, 2 or 3 means that x solves one of the problems ! Ax = b, min norm(Ax - b) or damped least squares. ! 4 means the limit on cond(A) was reached. ! 5 means the limit on iterations was reached. ! * 07 Dec 1994: Keep track of dxmax = max_k norm( phi_k * d_k ). ! So far, this is just printed at the end. ! A large value (relative to norm(x)) indicates ! significant cancellation in forming ! x = D*f = sum( phi_k * d_k ). ! A large column of D need NOT be serious if the ! corresponding phi_k is small. ! * 27 Dec 1994: Include estimate of alfa_opt in iteration log. ! alfa_opt is the optimal scale factor for the ! residual in the "augmented system", as described by ! A. Bjorck (1992), ! Pivoting and stability in the augmented system method, ! in D. F. Griffiths and G. A. Watson (eds.), ! "Numerical Analysis 1991", ! Proceedings of the 14th Dundee Conference, ! Pitman Research Notes in Mathematics 260, ! Longman Scientific and Technical, Harlow, Essex, 1992. ! * 12 Nov 2019 : Jacob Williams : significant refactoring into modern Fortran. ! !### Author ! * Michael A. Saunders, Dept of Operations Research, Stanford University ! !@note The number of iterations required by LSQR will usually decrease ! if the computation is performed in higher precision.

      • acheck

        public void acheck​(LSQR.lsqr_solver_ez me,
                   int m,
                   int n,
                   java.io.RandomAccessFile raFile,
                   double eps,
                   double[] v,
                   double[] w,
                   double[] x,
                   double[] y,
                   int[] inform)
        !> ! Checks the two modes of aprod for [[lsqr]]. ! ! acheck may be called to test the user-written subroutine ! aprod required by LSQR and CRAIG. For some m x n matrix A, ! aprod with mode = 1 and 2 supplies LSQR and CRAIG with products ! of the form ! y := y + Ax and x := x + A'y ! respectively, where A' means A(transpose). ! acheck tries to verify that A and A' refer to the same matrix. ! !### Method ! We cook up some "unlikely" vectors x and y of unit length ! and test if y'(y + Ax) = x'(x + A'y). ! !### History ! * 04 Sep 1991 Initial design and code. ! Michael Saunders, Dept of Operations Research, ! Stanford University ! * 10 Feb 1992 aprod and eps added as parameters. ! tol defined via power.

      • xcheck

        public void xcheck​(LSQR.lsqr_solver_ez me,
                   int m,
                   int n,
                   java.io.RandomAccessFile raFile,
                   double anorm,
                   double damp,
                   double eps,
                   double[] b,
                   double[] u,
                   double[] v,
                   double[] w,
                   double[] x,
                   int[] inform,
                   double[] test1,
                   double[] test2,
                   double[] test3)
        !> ! Tests if `x` solves a certain least-squares problem. ! ! xcheck computes residuals and norms associated with the ! vector x and the least-squares problem solved by LSQR or CRAIG. ! It determines whether x seems to be a solution to any of three ! possible systems: ! ! 1. Ax = b ! 2. min norm(Ax - b) ! 3. min norm(Ax - b)^2 + damp^2 * norm(x)^2 ! !### History ! * 07 Feb 1992 Initial design and code. ! Michael Saunders, Dept of Operations Research, ! Stanford University.

      • lsqrtest_ez

        public void lsqrtest_ez()

      • test_1

        public void test_1()

      • test_2

        public void test_2()

      • ddot

        double ddot​(int n,
            double[] dx,
            int incx,
            double[] dy,
            int incy)

      • d2norm

        double d2norm​(double a,
              double b)
        !> ! Returns \( \sqrt{ a^2 + b^2 } \) with precautions to avoid overflow. ! !### History ! * 21 Mar 1990: First version.

      • dnrm2

        double dnrm2​(int n,
             double[] x,
             int incx)

      • dscal

        int dscal​(int n,
          double da,
          double[] dx,
          int incx)

      • dcopy

        int dcopy​(int n,
          double[] dx,
          int incx,
          double[] dy,
          int incy)

      • lsqr__test

        public void lsqr__test()

      • test

        public void test​(LSQR.lsqr_solver_ez me,
                 int m,
                 int n,
                 int nduplc,
                 int npower,
                 double damp)
        !> ! This is an example driver routine for running LSQR. ! It generates a test problem, solves it, and examines the results. ! Note that subroutine aprod must be declared external ! if it is used only in the call to LSQR (and acheck). ! !### History ! * 1982---1991: Various versions implemented. ! * 04 Sep 1991: "wantse" added to argument list of LSQR, ! making standard errors optional. ! * 10 Feb 1992: Revised for use with lsqrchk fortran. ! * 31 Mar 2005: changed atol = eps**0.666667 to eps*0.99 ! to increase accuracy of the solution. LSQR appears to be ! successful on all 18 test problems except 5 and 6 ! (which are over-determined and too ill-conditioned to ! permit any correct digits). ! The output from an Intel Xeon system with g77 is in LSQR.LIS. ! The two "appears to have failed" messages are no cause for alarm. ! Michael Saunders, Dept of Operations Research, Stanford University.

      • lstp

        public void lstp​(LSQR.lsqr_solver_ez me,
                 int m,
                 int n,
                 int maxmn,
                 int minmn,
                 int nduplc,
                 int npower,
                 double damp,
                 double[] x,
                 double[] b,
                 double[] d,
                 double[] hy,
                 double[] hz,
                 double[] w,
                 double[] acond,
                 double[] rnorm)
        !> ! lstp generate a sparse least-squares test problem of the form ! ( A )*x = ( b ) ! ( damp*I ) ( 0 ) ! for solution by LSQR, or a sparse underdetermined system ! Ax + damp*s = b ! for solution by CRAIG. The matrix A is m by n and is ! constructed in the form A = HY*D*HZ, where D is an m by n ! diagonal matrix, and HY and HZ are Householder transformations. ! ! m and n may contain any positive values. ! The first 8 parameters are input to lstp. The last 8 are output. ! If m >= n or damp = 0, the true solution is x as given. ! Otherwise, x is modified to contain the true solution.

      • aprod1

        public void aprod1​(LSQR.lsqr_solver_ez me,
                   int m,
                   int n,
                   int maxmn,
                   int minmn,
                   double[] x,
                   double[] y,
                   double[] d,
                   double[] hy,
                   double[] hz,
                   double[] w)
        !> ! aprod1 computes y = y + A*x for subroutine aprod, ! where A is a test matrix of the form A = HY*D*HZ, ! and the latter matrices HY, D, HZ are represented by ! input vectors with the same name.

      • hprod

        public void hprod​(int n,
                  double[] hz,
                  double[] x,
                  double[] y)
        !> ! hprod applies a Householder transformation stored in hz ! to get y = ( I - 2*hz*hz(transpose) ) * x.