[Numpy-discussion] Problem in LinearAlgebra?

[Rob Hooft]

Nadav Horesh wrote:
> Many unstable problems have a stable solution if you choose the right 
> algorithm. The question is if somewhere the developers decided to switch 
> the equations solver, or there is a real bug. I hope that one of the 
> developers will reply in this forum. I will try to look at it also, 
> since it is a core component in the linear algebra package. Also if you 
> suspect that your work-around is useful --- please post it here!
> 

The workaround is to use "generalized_inverse" instead of 
"solve_linear_equations".

The changes in the latter routine since 17.1.2 are:

@@ -269,18 +408,37 @@
      bstar = Numeric.zeros((ldb,n_rhs),t)
      bstar[:b.shape[0],:n_rhs] = copy.copy(b)
      a,bstar = _castCopyAndTranspose(t, a, bstar)
-    lwork = 8*min(n,m) + max([2*min(m,n),max(m,n),n_rhs])
      s = Numeric.zeros((min(m,n),),real_t)
-    work = Numeric.zeros((lwork,), t)
+    nlvl = max( 0, int( math.log( float(min( m,n ))/2. ) ) + 1 )
+    iwork = Numeric.zeros((3*min(m,n)*nlvl+11*min(m,n),), 'l')
      if _array_kind[t] == 1: # Complex routines take different arguments
-        lapack_routine = lapack_lite.zgelss
-        rwork = Numeric.zeros((5*min(m,n)-1,), real_t)
+        lapack_routine = lapack_lite.zgelsd
+        lwork = 1
+        rwork = Numeric.zeros((lwork,), real_t)
+        work = Numeric.zeros((lwork,),t)
          results = lapack_routine( m, n, n_rhs, a, m, bstar,ldb , s, rcond,
-                        0,work,lwork,rwork,0 )
+                        0,work,-1,rwork,iwork,0 )
+        lwork = int(abs(work[0]))
+        rwork = Numeric.zeros((lwork,),real_t)
+        a_real = Numeric.zeros((m,n),real_t)
+        bstar_real = Numeric.zeros((ldb,n_rhs,),real_t)
+        results = lapack_lite.dgelsd( m, n, n_rhs, a_real, m, 
bstar_real,ldb , s, rcond,
+                        0,rwork,-1,iwork,0 )
+        lrwork = int(rwork[0])
+        work = Numeric.zeros((lwork,), t)
+        rwork = Numeric.zeros((lrwork,), real_t)
+        results = lapack_routine( m, n, n_rhs, a, m, bstar,ldb , s, rcond,
+                        0,work,lwork,rwork,iwork,0 )
      else:
-        lapack_routine = lapack_lite.dgelss
+        lapack_routine = lapack_lite.dgelsd
+        lwork = 1
+        work = Numeric.zeros((lwork,), t)
+        results = lapack_routine( m, n, n_rhs, a, m, bstar,ldb , s, rcond,
+                        0,work,-1,iwork,0 )
+        lwork = int(work[0])
+        work = Numeric.zeros((lwork,), t)
          results = lapack_routine( m, n, n_rhs, a, m, bstar,ldb , s, rcond,
-                        0,work,lwork,0 )
+                        0,work,lwork,iwork,0 )
      if results['info'] > 0:
          raise LinAlgError, 'SVD did not converge in Linear Least Squares'
      resids = Numeric.array([],t)

I'm not deep enough into this to know where the new version goes wrong.

Regards,

Rob Hooft

-- 
Rob W.W. Hooft  ||  rob at hooft.net  ||  http://www.hooft.net/people/rob/