<s> '''<STR_LIT>''' <EOL> from __future__ import print_function <EOL> from statsmodels . compat . python import lmap , range <EOL> import numpy as np <EOL> from scipy . linalg import svd as decomp_svd <EOL> from scipy . linalg . lapack import get_lapack_funcs <EOL> from numpy import asarray , zeros , sum , conjugate , dot , transpose <EOL> import numpy <EOL> from numpy import asarray_chkfinite , single <EOL> from numpy . linalg import LinAlgError <EOL> def lstsq ( a , b , cond = None , overwrite_a = <NUM_LIT:0> , overwrite_b = <NUM_LIT:0> ) : <EOL> """<STR_LIT>""" <EOL> a1 , b1 = lmap ( asarray_chkfinite , ( a , b ) ) <EOL> if a1 . ndim != <NUM_LIT:2> : <EOL> raise ValueError ( '<STR_LIT>' ) <EOL> m , n = a1 . shape <EOL> if b1 . ndim == <NUM_LIT:2> : <EOL> nrhs = b1 . shape [ <NUM_LIT:1> ] <EOL> else : <EOL> nrhs = <NUM_LIT:1> <EOL> if m != b1 . shape [ <NUM_LIT:0> ] : <EOL> raise ValueError ( '<STR_LIT>' ) <EOL> gelss , = get_lapack_funcs ( ( '<STR_LIT>' , ) , ( a1 , b1 ) ) <EOL> if n > m : <EOL> b2 = zeros ( ( n , nrhs ) , dtype = gelss . dtype ) <EOL> if b1 . ndim == <NUM_LIT:2> : <EOL> b2 [ : m , : ] = b1 <EOL> else : <EOL> b2 [ : m , <NUM_LIT:0> ] = b1 <EOL> b1 = b2 <EOL> overwrite_a = overwrite_a or ( a1 is not a and not hasattr ( a , '<STR_LIT>' ) ) <EOL> overwrite_b = overwrite_b or ( b1 is not b and not hasattr ( b , '<STR_LIT>' ) ) <EOL> if gelss . module_name [ : <NUM_LIT:7> ] == '<STR_LIT>' : <EOL> work = gelss ( a1 , b1 , lwork = - <NUM_LIT:1> ) [ <NUM_LIT:4> ] <EOL> lwork = work [ <NUM_LIT:0> ] . real . astype ( np . int ) <EOL> v , x , s , rank , work , info = gelss ( <EOL> a1 , b1 , cond = cond , lwork = lwork , overwrite_a = overwrite_a , <EOL> overwrite_b = overwrite_b ) <EOL> else : <EOL> raise NotImplementedError ( '<STR_LIT>' % <EOL> gelss . module_name ) <EOL> if info > <NUM_LIT:0> : <EOL> raise LinAlgError ( "<STR_LIT>" ) <EOL> if info < <NUM_LIT:0> : <EOL> raise ValueError ( '<STR_LIT>' <EOL> '<STR_LIT>' % - info ) <EOL> resids = asarray ( [ ] , dtype = x . dtype ) <EOL> if n < m : <EOL> x1 = x [ : n ] <EOL> if rank == n : <EOL> resids = sum ( x [ n : ] ** <NUM_LIT:2> , axis = <NUM_LIT:0> ) <EOL> x = x1 <EOL> return x , resids , rank , s <EOL> def pinv ( a , cond = None , rcond = None ) : <EOL> """<STR_LIT>""" <EOL> a = asarray_chkfinite ( a ) <EOL> b = numpy . identity ( a . shape [ <NUM_LIT:0> ] , dtype = a . dtype ) <EOL> if rcond is not None : <EOL> cond = rcond <EOL> return lstsq ( a , b , cond = cond ) [ <NUM_LIT:0> ] <EOL> eps = numpy . finfo ( float ) . eps <EOL> feps = numpy . finfo ( single ) . eps <EOL> _array_precision = { '<STR_LIT:f>' : <NUM_LIT:0> , '<STR_LIT:d>' : <NUM_LIT:1> , '<STR_LIT:F>' : <NUM_LIT:0> , '<STR_LIT:D>' : <NUM_LIT:1> } <EOL> def pinv2 ( a , cond = None , rcond = None ) : <EOL> """<STR_LIT>""" <EOL> a = asarray_chkfinite ( a ) <EOL> u , s , vh = decomp_svd ( a ) <EOL> t = u . dtype . char <EOL> if rcond is not None : <EOL> cond = rcond <EOL> if cond in [ None , - <NUM_LIT:1> ] : <EOL> cond = { <NUM_LIT:0> : feps * <NUM_LIT> , <NUM_LIT:1> : eps * <NUM_LIT> } [ _array_precision [ t ] ] <EOL> m , n = a . shape <EOL> cutoff = cond * numpy . maximum . reduce ( s ) <EOL> psigma = zeros ( ( m , n ) , t ) <EOL> for i in range ( len ( s ) ) : <EOL> if s [ i ] > cutoff : <EOL> psigma [ i , i ] = <NUM_LIT:1.0> / conjugate ( s [ i ] ) <EOL> return transpose ( conjugate ( dot ( dot ( u , psigma ) , vh ) ) ) <EOL> def logdet_symm ( m , check_symm = False ) : <EOL> """<STR_LIT>""" <EOL> from scipy import linalg <EOL> if check_symm : <EOL> if not np . all ( m == m . T ) : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> c , _ = linalg . cho_factor ( m , lower = True ) <EOL> return <NUM_LIT:2> * np . sum ( np . log ( c . diagonal ( ) ) ) <EOL> def stationary_solve ( r , b ) : <EOL> """<STR_LIT>""" <EOL> db = r [ <NUM_LIT:0> : <NUM_LIT:1> ] <EOL> dim = b . ndim <EOL> if b . ndim == <NUM_LIT:1> : <EOL> b = b [ : , None ] <EOL> x = b [ <NUM_LIT:0> : <NUM_LIT:1> , : ] <EOL> for j in range ( <NUM_LIT:1> , len ( b ) ) : <EOL> rf = r [ <NUM_LIT:0> : j ] [ : : - <NUM_LIT:1> ] <EOL> a = ( b [ j , : ] - np . dot ( rf , x ) ) / ( <NUM_LIT:1> - np . dot ( rf , db [ : : - <NUM_LIT:1> ] ) ) <EOL> z = x - np . outer ( db [ : : - <NUM_LIT:1> ] , a ) <EOL> x = np . concatenate ( ( z , a [ None , : ] ) , axis = <NUM_LIT:0> ) <EOL> if j == len ( b ) - <NUM_LIT:1> : <EOL> break <EOL> rn = r [ j ] <EOL> a = ( rn - np . dot ( rf , db ) ) / ( <NUM_LIT:1> - np . dot ( rf , db [ : : - <NUM_LIT:1> ] ) ) <EOL> z = db - a * db [ : : - <NUM_LIT:1> ] <EOL> db = np . concatenate ( ( z , np . r_ [ a ] ) ) <EOL> if dim == <NUM_LIT:1> : <EOL> x = x [ : , <NUM_LIT:0> ] <EOL> return x <EOL> if __name__ == '<STR_LIT:__main__>' : <EOL> a0 = np . random . randn ( <NUM_LIT:100> , <NUM_LIT:10> ) <EOL> b0 = a0 . sum ( <NUM_LIT:1> ) [ : , None ] + np . random . randn ( <NUM_LIT:100> , <NUM_LIT:3> ) <EOL> lstsq ( a0 , b0 ) <EOL> pinv ( a0 ) <EOL> pinv2 ( a0 ) <EOL> x = pinv ( a0 ) <EOL> x2 = scipy . linalg . pinv ( a0 ) <EOL> print ( np . max ( np . abs ( x - x2 ) ) ) <EOL> x = pinv2 ( a0 ) <EOL> x2 = scipy . linalg . pinv2 ( a0 ) <EOL> print ( np . max ( np . abs ( x - x2 ) ) ) </s>
