[PYTHON MATRIX-SIG] Linear Least Squares for LinAlg.py

[Doug Heisterkamp]

Hi,

I've written a linear least square function to add to LinAlg.py.  Let
me know if you find any problems with it.  

Doug Heisterkamp
drh@cherokee.unl.edu



# Linear Least Squares 
	
def solveLinearLeastSquares(a, b, rcond=1.e-10):
    """solveLinearLeastSquares(a,b) returns x,resids,rank,s  
where x minimizes 2-norm(|b - Ax|)
      resids is the sum square residuals
      rank is the rank of A
      s is an rank of the singual values of A in desending order
If b is a matrix then x is also a matrix with corresponding columns.
If the rank of A is less than the number of columns of A or greater than
the numer of rows, then residuals will be returned as an empty array
otherwise resids = sum((b-dot(A,x)**2).
Singular values less than s[0]*rcond are treated as zero.
"""
    one_eq = len(b.shape) == 1
    if one_eq:
	b = b[:, Numeric.NewAxis]
    _assertRank2(a, b)
#   _assertSquareness(a)
    m  = a.shape[0]
    n  = a.shape[1]
    n_rhs = b.shape[1]
    ldb = max(n,m)
    if m != b.shape[0]:
        raise LinAlgError, 'Incompatible dimensions'
    t =_commonType(a, b)
    real_t = _array_type[0][_array_precision[t]]
    lapack_routine = _findLapackRoutine('gelss', t)
    bstar = Numeric.zeros((ldb,n_rhs),t)
    bstar[:b.shape[0],:n_rhs] = cp(b)
    a,bstar = _castCopyAndTranspose(t, a, bstar)
#   minimum value of lwork:3*min(n,m) + max([2*min(m,n),max(m,n),n_rhs])
    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)
#   Note: in next version rcond will be passed directly to lapack module
#         by for now wrap it in an array
    arcond = Numeric.array([rcond],real_t)
    if _array_kind[t] == 1: # Complex routines take different arguments
	rwork = Numeric.zeros((5*min(m,n)-1,), real_t)
        results = lapack_routine( m, n, n_rhs, a, m, bstar,ldb , s, arcond,
			0,work,lwork,rwork,0 )
    else:
        results = lapack_routine( m, n, n_rhs, a, m, bstar,ldb , s, arcond,
			0,work,lwork,0 )
    if results['info'] > 0:
	raise LinAlgError, 'SVD did not converge in Linear Least Squares'
    resids = Numeric.array([],t)
    if one_eq:
	x = cp(Numeric.ravel(bstar)[:n])
        if (results['rank']==n) and (m>n):
            resids = Numeric.array([Numeric.sum((Numeric.ravel(bstar)[n:])**2)])
    else:
	x = cp(transpose(bstar)[:n,:])
        if (results['rank']==n) and (m>n):
            resids = cp(Numeric.sum((transpose(bstar)[n:,:])**2))
    return x,resids,results['rank'],cp(s[:min(n,m)]) 


=================
MATRIX-SIG  - SIG on Matrix Math for Python

send messages to: matrix-sig@python.org
administrivia to: matrix-sig-request@python.org
=================