c464e5d30a46bb15587c50b39ce9ab9cafec7e28
[stefan/supreme.git] / supreme / resolve / lsqr.py
1 """Sparse Equations and Least Squares.
2
3 The original Fortran code was written by C. C. Paige and M. A. Saunders as
4 described in
5
6 C. C. Paige and M. A. Saunders, LSQR: An algorithm for sparse linear
7 equations and sparse least squares, TOMS 8(1), 43--71 (1982).
8
9 C. C. Paige and M. A. Saunders, Algorithm 583; LSQR: Sparse linear
10 equations and least-squares problems, TOMS 8(2), 195--209 (1982).
11
12 It is licensed under the following BSD license:
13
14 Copyright (c) 2006, Systems Optimization Laboratory
15 All rights reserved.
16
17 Redistribution and use in source and binary forms, with or without
18 modification, are permitted provided that the following conditions are
19 met:
20
21     * Redistributions of source code must retain the above copyright
22       notice, this list of conditions and the following disclaimer.
23
24     * Redistributions in binary form must reproduce the above
25       copyright notice, this list of conditions and the following
26       disclaimer in the documentation and/or other materials provided
27       with the distribution.
28
29     * Neither the name of Stanford University nor the names of its
30       contributors may be used to endorse or promote products derived
31       from this software without specific prior written permission.
32
33 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
34 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
35 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
36 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
37 OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
38 SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
39 LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
40 DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
41 THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
42 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
43 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
44
45 The Fortran code was translated to Python for use in CVXOPT by Jeffery
46 Kline with contributions by Mridul Aanjaneya and Bob Myhill.
47
48 Adapted for SciPy by Stefan van der Walt.
49
50 """
51
52 __all__ = ['lsqr']
53
54 import numpy as np
55 from math import sqrt
56 from scipy.sparse.linalg import aslinearoperator
57
58 def _sym_ortho(a,b):
59     """
60     Jeffery Kline noted: I added the routine 'SymOrtho' for numerical
61     stability. This is recommended by S.-C. Choi in [1]_. It removes
62     the unpleasant potential of ``1/eps`` in some important places
63     (see, for example text following "Compute the next
64     plane rotation Qk" in minres_py).
65
66     References
67     ----------
68     .. [1] S.-C. Choi, "Iterative Methods for Singular Linear Equations
69            and Least-Squares Problems", Dissertation,
70            http://www.stanford.edu/group/SOL/dissertations/sou-cheng-choi-thesis.pdf
71
72     """
73     aa = abs(a)
74     ab = abs(b)
75     if b == 0.:
76         s = 0.
77         r = aa
78         if aa == 0.:
79             c = 1.
80         else:
81             c = a/aa
82     elif a == 0.:
83         c = 0.
84         s = b / ab
85         r = ab
86     elif ab >= aa:
87         sb = 1
88         if b < 0: sb=-1
89         tau = a/b
90         s = sb * (1 + tau**2)**-0.5
91         c = s * tau
92         r = b / s
93     elif aa > ab:
94         sa = 1
95         if a < 0: sa = -1
96         tau = b / a
97         c = sa * (1 + tau**2)**-0.5
98         s = c * tau
99         r = a / c
100
101     return c, s, r
102
103
104 def lsqr(A, b, damp=0.0, atol=1e-8, btol=1e-8, conlim=1e8,
105          iter_lim=None, show=False, calc_var=False):
106     """Find the least-squares solution to a large, sparse, linear system
107     of equations.
108
109     The function solves ``Ax = b``  or  ``min ||b - Ax||^2`` or
110     ``min ||Ax - b||^2 + d^2 ||x||^2.
111
112     The matrix A may be square or rectangular (over-determined or
113     under-determined), and may have any rank.
114
115     ::
116
117       1. Unsymmetric equations --    solve  A*x = b
118
119       2. Linear least squares  --    solve  A*x = b
120                                      in the least-squares sense
121
122       3. Damped least squares  --    solve  (   A    )*x = ( b )
123                                             ( damp*I )     ( 0 )
124                                      in the least-squares sense
125
126     Parameters
127     ----------
128     A : LinearOperator or equivalent
129         A representation of an mxn matrix.  It is required that
130         the linear operator can produce ``Ax`` and ``A.T x``.
131     b : (m,) ndarray
132         Right-hand side vector ``b``.
133     damp : float
134         Damping coefficient.
135     atol, btol : float
136         Stopping tolerances. If both are 1.0e-9 (say), the final
137         residual norm should be accurate to about 9 digits.  (The
138         final x will usually have fewer correct digits, depending on
139         cond(A) and the size of damp.)
140     conlim : float
141         Another stopping tolerance.  lsqr terminates if an estimate of
142         ``cond(A)`` exceeds `conlim`.  For compatible systems ``Ax =
143         b``, `conlim` could be as large as 1.0e+12 (say).  For
144         least-squares problems, conlim should be less than 1.0e+8.
145         Maximum precision can be obtained by setting ``atol = btol =
146         conlim = zero``, but the number of iterations may then be
147         excessive.
148     iter_lim : int
149         Explicit limitation on number of iterations (for safety).
150     show : bool
151         Display an iteration log.
152     calc_var : bool
153         Whether to estimate diagonals of ``(A'A + damp^2*I)^{-1}``.
154
155     Returns
156     -------
157     x : ndarray of float
158         The final solution.
159     istop : int
160         Gives the reason for termination.
161         1 means x is an approximate solution to Ax = b.
162         2 means x approximately solves the least-squares problem.
163     itn : int
164         Iteration number upon termination.
165     r1norm : float
166         ``norm(r)``, where ``r = b - Ax``.
167     r2norm : float
168         ``sqrt( norm(r)^2  +  damp^2 * norm(x)^2 )``.  Equal to `r1norm` if
169         ``damp == 0``.
170     anorm : float
171         Estimate of Frobenius norm of ``Abar = [[A]; [damp*I]]``.
172     acond : float
173         Estimate of ``cond(Abar)``.
174     arnorm : float
175         Estimate of ``norm(A'*r - damp^2*x)``.
176     xnorm : float
177         ``norm(x)``
178     var : ndarray of float
179         If ``calc_var`` is True, estimates all diagonals of
180         ``(A'A)^{-1}`` (if ``damp == 0``) or more generally ``(A'A +
181         damp^2*I)^{-1}``.  This is well defined if A has full column
182         rank or ``damp > 0``.  (Not sure what var means if ``rank(A)
183         < n`` and ``damp = 0.``)
184
185     Notes
186     -----
187     LSQR uses an iterative method to approximate the solution.  The
188     number of iterations required to reach a certain accuracy depends
189     strongly on the scaling of the problem.  Poor scaling of the rows
190     or columns of A should therefore be avoided where possible.
191
192     For example, in problem 1 the solution is unaltered by
193     row-scaling.  If a row of A is very small or large compared to
194     the other rows of A, the corresponding row of ( A  b ) should be
195     scaled up or down.
196
197     In problems 1 and 2, the solution x is easily recovered
198     following column-scaling.  Unless better information is known,
199     the nonzero columns of A should be scaled so that they all have
200     the same Euclidean norm (e.g., 1.0).
201
202     In problem 3, there is no freedom to re-scale if damp is
203     nonzero.  However, the value of damp should be assigned only
204     after attention has been paid to the scaling of A.
205
206     The parameter damp is intended to help regularize
207     ill-conditioned systems, by preventing the true solution from
208     being very large.  Another aid to regularization is provided by
209     the parameter acond, which may be used to terminate iterations
210     before the computed solution becomes very large.
211
212     If some initial estimate ``x0`` is known and if ``damp == 0``,
213     one could proceed as follows:
214
215       1. Compute a residual vector ``r0 = b - A*x0``.
216       2. Use LSQR to solve the system  ``A*dx = r0``.
217       3. Add the correction dx to obtain a final solution ``x = x0 + dx``.
218
219     This requires that ``x0`` be available before and after the call
220     to LSQR.  To judge the benefits, suppose LSQR takes k1 iterations
221     to solve A*x = b and k2 iterations to solve A*dx = r0.
222     If x0 is "good", norm(r0) will be smaller than norm(b).
223     If the same stopping tolerances atol and btol are used for each
224     system, k1 and k2 will be similar, but the final solution x0 + dx
225     should be more accurate.  The only way to reduce the total work
226     is to use a larger stopping tolerance for the second system.
227     If some value btol is suitable for A*x = b, the larger value
228     btol*norm(b)/norm(r0)  should be suitable for A*dx = r0.
229
230     Preconditioning is another way to reduce the number of iterations.
231     If it is possible to solve a related system ``M*x = b``
232     efficiently, where M approximates A in some helpful way (e.g. M -
233     A has low rank or its elements are small relative to those of A),
234     LSQR may converge more rapidly on the system ``A*M(inverse)*z =
235     b``, after which x can be recovered by solving M*x = z.
236
237     If A is symmetric, LSQR should not be used!
238
239     Alternatives are the symmetric conjugate-gradient method (cg)
240     and/or SYMMLQ.  SYMMLQ is an implementation of symmetric cg that
241     applies to any symmetric A and will converge more rapidly than
242     LSQR.  If A is positive definite, there are other implementations
243     of symmetric cg that require slightly less work per iteration than
244     SYMMLQ (but will take the same number of iterations).
245
246     References
247     ----------
248     .. [1] C. C. Paige and M. A. Saunders (1982a).
249            "LSQR: An algorithm for sparse linear equations and
250            sparse least squares", ACM TOMS 8(1), 43-71.
251     .. [2] C. C. Paige and M. A. Saunders (1982b).
252            "Algorithm 583.  LSQR: Sparse linear equations and least
253            squares problems", ACM TOMS 8(2), 195-209.
254     .. [3] M. A. Saunders (1995).  "Solution of sparse rectangular
255            systems using LSQR and CRAIG", BIT 35, 588-604.
256
257     """
258     A = aslinearoperator(A)
259     b = b.squeeze()
260
261     m, n = A.shape
262     if iter_lim is None: iter_lim = 2 * n
263     var = np.zeros(n)
264
265     msg=('The exact solution is  x = 0                              ',
266          'Ax - b is small enough, given atol, btol                  ',
267          'The least-squares solution is good enough, given atol     ',
268          'The estimate of cond(Abar) has exceeded conlim            ',
269          'Ax - b is small enough for this machine                   ',
270          'The least-squares solution is good enough for this machine',
271          'Cond(Abar) seems to be too large for this machine         ',
272          'The iteration limit has been reached                      ');
273
274     if show:
275         print ' '
276         print 'LSQR            Least-squares solution of  Ax = b'
277         str1 = 'The matrix A has %8g rows  and %8g cols' % (m, n)
278         str2 = 'damp = %20.14e   calc_var = %8g' % (damp, calc_var)
279         str3 = 'atol = %8.2e                 conlim = %8.2e'%( atol, conlim)
280         str4 = 'btol = %8.2e               iter_lim = %8g'  %( btol, iter_lim)
281         print str1
282         print str2
283         print str3
284         print str4
285
286     itn = 0
287     istop = 0
288     nstop = 0
289     ctol = 0
290     if conlim > 0: ctol = 1/conlim
291     anorm = 0
292     acond = 0
293     dampsq = damp**2
294     ddnorm = 0
295     res2 = 0
296     xnorm = 0
297     xxnorm = 0
298     z = 0
299     cs2 = -1
300     sn2 = 0
301
302     """
303     Set up the first vectors u and v for the bidiagonalization.
304     These satisfy  beta*u = b,  alfa*v = A'u.
305     """
306     __xm = np.zeros(m) # a matrix for temporary holding
307     __xn = np.zeros(n) # a matrix for temporary holding
308     v = np.zeros(n)
309     u = b
310     x = np.zeros(n)
311     alfa = 0
312     beta = np.linalg.norm(u)
313     w = np.zeros(n)
314
315     if beta > 0:
316         u = (1/beta) * u
317         v = A.rmatvec(u)
318         alfa = np.linalg.norm(v)
319
320     if alfa > 0:
321         v = (1/alfa) * v
322         w = v.copy()
323
324     rhobar = alfa
325     phibar = beta
326     bnorm = beta
327     rnorm = beta
328     r1norm = rnorm
329     r2norm = rnorm
330
331     # Reverse the order here from the original matlab code because
332     # there was an error on return when arnorm==0
333     arnorm = alfa * beta
334     if arnorm == 0:
335         print msg[0];
336         return x, istop, itn, r1norm, r2norm, anorm, acond, arnorm, xnorm, var
337
338     head1  = '   Itn      x[0]       r1norm     r2norm ';
339     head2  = ' Compatible    LS      Norm A   Cond A';
340
341     if show:
342         print ' '
343         print head1, head2
344         test1  = 1;             test2  = alfa / beta;
345         str1   = '%6g %12.5e'    %(    itn,   x[0] );
346         str2   = ' %10.3e %10.3e'%( r1norm, r2norm );
347         str3   = '  %8.1e %8.1e' %(  test1,  test2 );
348         print str1, str2, str3
349
350     # Main iteration loop.
351     while itn < iter_lim:
352         itn = itn + 1
353         """
354         %     Perform the next step of the bidiagonalization to obtain the
355         %     next  beta, u, alfa, v.  These satisfy the relations
356         %                beta*u  =  a*v   -  alfa*u,
357         %                alfa*v  =  A'*u  -  beta*v.
358         """
359         u = A.matvec(v) - alfa * u
360         beta = np.linalg.norm(u)
361
362         if beta > 0:
363             u = (1/beta) * u
364             anorm = sqrt(anorm**2 + alfa**2 + beta**2 + damp**2)
365             v = A.rmatvec(u) - beta * v
366             alfa  = np.linalg.norm(v)
367             if alfa > 0:
368                 v = (1 / alfa) * v
369
370         # Use a plane rotation to eliminate the damping parameter.
371         # This alters the diagonal (rhobar) of the lower-bidiagonal matrix.
372         rhobar1 = sqrt(rhobar**2 + damp**2)
373         cs1 = rhobar / rhobar1
374         sn1 = damp / rhobar1
375         psi = sn1 * phibar
376         phibar = cs1 * phibar
377
378         # Use a plane rotation to eliminate the subdiagonal element (beta)
379         # of the lower-bidiagonal matrix, giving an upper-bidiagonal matrix.
380         cs, sn, rho = _sym_ortho(rhobar1, beta)
381
382         theta = sn * alfa
383         rhobar = -cs * alfa
384         phi = cs * phibar
385         phibar = sn * phibar
386         tau = sn * phi
387
388         # Update x and w.
389         t1 = phi / rho
390         t2 = -theta / rho
391         dk = (1 / rho) * w
392
393         x = x + t1 * w
394         w = v + t2 * w
395         ddnorm = ddnorm + np.linalg.norm(dk)**2
396
397         if calc_var:
398             var = var + dk**2
399
400         # Use a plane rotation on the right to eliminate the
401         # super-diagonal element (theta) of the upper-bidiagonal matrix.
402         # Then use the result to estimate norm(x).
403         delta = sn2 * rho
404         gambar = -cs2 * rho
405         rhs = phi - delta * z
406         zbar = rhs / gambar
407         xnorm = sqrt(xxnorm + zbar**2)
408         gamma = sqrt(gambar**2 +theta**2)
409         cs2 = gambar / gamma
410         sn2 = theta  / gamma
411         z = rhs / gamma
412         xxnorm = xxnorm  +  z**2
413
414         # Test for convergence.
415         # First, estimate the condition of the matrix  Abar,
416         # and the norms of  rbar  and  Abar'rbar.
417         acond = anorm * sqrt(ddnorm)
418         res1 = phibar**2
419         res2 = res2 + psi**2
420         rnorm = sqrt(res1 + res2)
421         arnorm = alfa * abs(tau)
422
423         # Distinguish between
424         #    r1norm = ||b - Ax|| and
425         #    r2norm = rnorm in current code
426         #           = sqrt(r1norm^2 + damp^2*||x||^2).
427         #    Estimate r1norm from
428         #    r1norm = sqrt(r2norm^2 - damp^2*||x||^2).
429         # Although there is cancellation, it might be accurate enough.
430         r1sq = rnorm**2 - dampsq * xxnorm
431         r1norm = sqrt(abs(r1sq))
432         if r1sq < 0:
433             r1norm = -r1norm
434         r2norm = rnorm
435
436         # Now use these norms to estimate certain other quantities,
437         # some of which will be small near a solution.
438         test1 = rnorm / bnorm
439         test2 = arnorm / (anorm * rnorm)
440         test3 = 1 / acond
441         t1 = test1 / (1 + anorm * xnorm / bnorm)
442         rtol = btol + atol *  anorm * xnorm / bnorm
443
444         # The following tests guard against extremely small values of
445         # atol, btol  or  ctol.  (The user may have set any or all of
446         # the parameters  atol, btol, conlim  to 0.)
447         # The effect is equivalent to the normal tests using
448         # atol = eps,  btol = eps,  conlim = 1/eps.
449         if itn >= iter_lim: istop = 7
450         if 1 + test3 <= 1: istop = 6
451         if 1 + test2 <= 1: istop = 5
452         if 1 + t1 <= 1: istop = 4
453
454         # Allow for tolerances set by the user.
455         if test3 <= ctol: istop = 3
456         if test2 <= atol: istop = 2
457         if test1 <= rtol: istop = 1
458
459         # See if it is time to print something.
460         prnt = False;
461         if n <= 40: prnt = True
462         if itn <= 10: prnt = True
463         if itn >= iter_lim-10: prnt = True
464         # if itn%10 == 0: prnt = True
465         if test3 <= 2*ctol: prnt = True
466         if test2 <= 10*atol: prnt = True
467         if test1 <= 10*rtol: prnt = True
468         if istop != 0: prnt = True
469
470         if prnt:
471             if show:
472                 str1 = '%6g %12.5e' % (itn, x[0])
473                 str2 = ' %10.3e %10.3e' % (r1norm, r2norm)
474                 str3 = '  %8.1e %8.1e' % (test1, test2)
475                 str4 = ' %8.1e %8.1e' % (anorm, acond)
476                 print str1, str2, str3, str4
477
478         if istop != 0: break
479
480     # End of iteration loop.
481     # Print the stopping condition.
482     if show:
483         print ' '
484         print 'LSQR finished'
485         print msg[istop]
486         print ' '
487         str1 = 'istop =%8g   r1norm =%8.1e' % (istop, r1norm)
488         str2 = 'anorm =%8.1e   arnorm =%8.1e' % (anorm, arnorm)
489         str3 = 'itn   =%8g   r2norm =%8.1e' % (itn, r2norm)
490         str4 = 'acond =%8.1e   xnorm  =%8.1e' % (acond, xnorm)
491         print str1+ '   ' + str2
492         print str3+ '   ' + str4
493         print ' '
494
495     return x, istop, itn, r1norm, r2norm, anorm, acond, arnorm, xnorm, var