BDS 761: Data Science and Machine Learning I
Topic 3: Core Libraries
This topic:¶
- BLAS
- Overview of major python packages
- Numpy practice
Readings:
- Harris, Charles R., et al. "Array programming with NumPy." nature 585.7825 (2020): 357-362.
- some python and linear algebra review material: https://www.keithdillon.com/index.php/bootcamp/
- netlib.org BLAS FAQ https://netlib.org/blas/faq.html
What is "Classical" Linear Algebra software?¶
- A term I just invented
- Originally single CPU, RAM-limited implementations of well-known algorithms (e.g. Golub & Van Loan "Matrix Computations")
- ...and higher-level packages based on these algorithms, inheriting the limitations
- Large-scale (more than 10k variable or so) only possible via sparse linear algebra algorithms
- Leading parallel tool is SIMD (single instruction multiple data) via math co-processors.
- More parallelism at top-level rather than bottom ("workers" running seperate processes on separate CPU's)
- Updates based on producing variants rather than complete replacement
1. Hardware Level (Machine Code / Registers)¶
Manually multiply and sum numbers stored in registers
; Pseudocode: compute 2*3 + 4*5
MOV R1, 2
MOV R2, 3
MUL R3, R1, R2 ; R3 = 6
MOV R4, 4
MOV R5, 5
MUL R6, R4, R5 ; R6 = 20
ADD R7, R3, R6 ; R7 = 26
In your computer memory it actually looks like this. Left is the binary data there, right is translation. This is an 8-bit computer.
00010001 00000010 ; MOV R1, 2
00010010 00000011 ; MOV R2, 3
00100011 00010010 ; MUL R3, R1, R2
00010100 00000100 ; MOV R4, 4
00010101 00000101 ; MOV R5, 5
00100110 01000101 ; MUL R6, R4, R5
00110111 00110110 ; ADD R7, R3, R6
3. C Code (Low-level Language)¶
Still close to hardware, but (mostly) portable across machines
int dot_product(int a[], int b[], int n) {
int sum = 0;
for (int i = 0; i < n; i++) {
sum += a[i] * b[i];
}
return sum;
}
Compilers convert this into assembly.
4. Linear Algebra Library (e.g., BLAS)¶
Wraps optimized native implementations
#include <cblas.h>
double result = cblas_ddot(n, a, 1, b, 1); // a·b
5. Machine Learning Framework (e.g., PyTorch)¶
Automatically optimized, abstracted away
import torch
a = torch.tensor([2.0, 4.0])
b = torch.tensor([3.0, 5.0])
result = torch.dot(a, b) # 2*3 + 4*5 = 6 + 20 = 26
Under the hood: this may call architecture-optimized code (AVX, GPU, etc.)
Summary of Abstraction Layers¶
| Level | Example |
|---|---|
| Hardware | CPU registers, instructions |
| Assembly | MOV, MUL |
| C code | loops, arrays |
| Linear algebra lib | cblas_ddot |
| ML Framework | torch.dot() |
Each layer hides the complexity of the one below it, allowing developers to focus on higher-level tasks without reimplementing basic operations.
I. BLAS+¶
BLAS = Basic Linear Algebra Subprograms¶
"Building block" algorithms for matrix algebra
- scaling, multiplying. "Level 1,2,3 functions"
- single and double precision, real and complex numbers
- generally don't take advantage of multiple threads or CPU's by chip maker
- Are often optimized for complex instructions however, e.g. SIMD
- Originally FORTRAN - one-based, columns first, like matlab
- various implementations, including parallelizable versions, C/C++ CUDA,
- http://www.netlib.org/blas/
- https://numpy.org/doc/stable/building/blas_lapack.html
More Implementations¶
- cuBLAS - NVIDIA https://developer.nvidia.com/cublas
- rocBLAS - AMD https://github.com/ROCm/rocBLAS?tab=readme-ov-file
- clBLAS, CLBlast - OpenCL https://github.com/CNugteren/CLBlast
- hipBLAS - a wrapper around other BLAS implementations https://github.com/ROCm/hipBLAS
- gotoBLAS, BLIS, OpenBLAS - optimized implementations for CPU https://en.wikipedia.org/wiki/GotoBLAS, https://www.openblas.net/, https://github.com/flame/blis
cuBLAS: Basic Linear Algebra on NVIDIA GPUs¶
https://developer.nvidia.com/cublas
https://stackoverflow.com/questions/25836003/normal-cuda-vs-cublas
BLAS Levels 1,2,3 ~ $O(n), O(n^2), O(n^3)$¶
In [24]:
import scipy.linalg.blas as blas
printcols(dir(blas),6)
HAS_ILP64 chemm daxpy dtrmv sspmv zhbmv __all__ chemv dcopy dtrsm sspr zhemm __builtins__ cher ddot dtrsv sspr2 zhemv __cached__ cher2 dgbmv dzasum sswap zher __doc__ cher2k dgemm dznrm2 ssymm zher2 __file__ cherk dgemv find_best_blas_type ssymv zher2k __loader__ chpmv dger functools ssyr zherk __name__ chpr dnrm2 get_blas_funcs ssyr2 zhpmv __package__ chpr2 drot icamax ssyr2k zhpr __spec__ crotg drotg idamax ssyrk zhpr2 _blas_alias cscal drotm isamax stbmv zrotg _cblas cspmv drotmg izamax stbsv zscal _fblas cspr dsbmv sasum stpmv zspmv _fblas_64 csrot dscal saxpy stpsv zspr _get_funcs csscal dspmv scasum strmm zswap _memoize_get_funcs cswap dspr scnrm2 strmv zsymm _np csymm dspr2 scopy strsm zsyr _type_conv csyr dswap sdot strsv zsyr2k _type_score csyr2k dsymm sgbmv zaxpy zsyrk caxpy csyrk dsymv sgemm zcopy ztbmv ccopy ctbmv dsyr sgemv zdotc ztbsv cdotc ctbsv dsyr2 sger zdotu ztpmv cdotu ctpmv dsyr2k snrm2 zdrot ztpsv cgbmv ctpsv dsyrk srot zdscal ztrmm cgemm ctrmm dtbmv srotg zgbmv ztrmv cgemv ctrmv dtbsv srotm zgemm ztrsm cgerc ctrsm dtpmv srotmg zgemv ztrsv cgeru ctrsv dtpsv ssbmv zgerc chbmv dasum dtrmm sscal zgeru
LAPACK = Linear Algebra Package¶
- Fortran, Uses BLAS as backend (can choose different implementations)
- http://www.netlib.org/lapack/
- https://github.com/Reference-LAPACK/lapack
- systems of simultaneous linear equations
- least-squares solutions of linear systems of equations
- eigenvalue problems
- singular value problems.
- associated matrix factorizations (LU, Cholesky, QR, SVD, Schur, generalized Schur)
- related computations such as reordering of the Schur factorizations and estimating condition numbers.
- Dense and banded matrices are handled, but not general sparse matrices.
- real and complex matrices
- both single and double precision.
- Matlab, and the Python tools we will use were essentially built on this
In [28]:
import scipy.linalg.lapack as lapack
printcols(dir(lapack),6)
HAS_ILP64 clartg dgetc2 dtrsyl ssbev zgttrs _DeprecatedImport claswp dgetrf dtrtri ssbevd zhbevd __all__ clauum dgetri dtrtrs ssbevx zhbevx __builtins__ cpbsv dgetri_lwork dtrttf ssfrk zhecon __cached__ cpbtrf dgetrs dtrttp sstebz zheequb __doc__ cpbtrs dgges dtzrzf sstein zheev __file__ cpftrf dggev dtzrzf_lwork sstemr zheev_lwork __loader__ cpftri dgglse flapack sstemr_lwork zheevd __name__ cpftrs dgglse_lwork get_lapack_funcs ssterf zheevd_lwork __package__ cpocon dgtsv ilaver sstev zheevr __spec__ cposv dgtsvx routine ssycon zheevr_lwork _check_work_float cposvx dgttrf sgbsv ssyconv zheevx _clapack cpotrf dgttrs sgbtrf ssyequb zheevx_lwork _compute_lwork cpotri dlamch sgbtrs ssyev zhegst _flapack cpotrs dlange sgebal ssyev_lwork zhegv _flapack_64 cppcon dlarf sgecon ssyevd zhegv_lwork _get_funcs cppsv dlarfg sgeequ ssyevd_lwork zhegvd _int32_max cpptrf dlartg sgeequb ssyevr zhegvx _int64_max cpptri dlasd4 sgees ssyevr_lwork zhegvx_lwork _lapack_alias cpptrs dlaswp sgeev ssyevx zhesv _memoize_get_funcs cpstf2 dlauum sgeev_lwork ssyevx_lwork zhesv_lwork _np cpstrf dorcsd sgehrd ssygst zhesvx cgbsv cpteqr dorcsd_lwork sgehrd_lwork ssygv zhesvx_lwork cgbtrf cptsv dorghr sgejsv ssygv_lwork zhetrd cgbtrs cptsvx dorghr_lwork sgels ssygvd zhetrd_lwork cgebal cpttrf dorgqr sgels_lwork ssygvx zhetrf cgecon cpttrs dorgrq sgelsd ssygvx_lwork zhetrf_lwork cgeequ crot dormqr sgelsd_lwork ssysv zhfrk cgeequb csycon dormrz sgelss ssysv_lwork zlange cgees csyconv dormrz_lwork sgelss_lwork ssysvx zlarf cgeev csyequb dpbsv sgelsy ssysvx_lwork zlarfg cgeev_lwork csysv dpbtrf sgelsy_lwork ssytf2 zlartg cgehrd csysv_lwork dpbtrs sgemqrt ssytrd zlaswp cgehrd_lwork csysvx dpftrf sgeqp3 ssytrd_lwork zlauum cgels csysvx_lwork dpftri sgeqrf ssytrf zpbsv cgels_lwork csytf2 dpftrs sgeqrf_lwork ssytrf_lwork zpbtrf cgelsd csytrf dpocon sgeqrfp stbtrs zpbtrs cgelsd_lwork csytrf_lwork dposv sgeqrfp_lwork stfsm zpftrf cgelss ctbtrs dposvx sgeqrt stfttp zpftri cgelss_lwork ctfsm dpotrf sgerqf stfttr zpftrs cgelsy ctfttp dpotri sgesc2 stgexc zpocon cgelsy_lwork ctfttr dpotrs sgesdd stgsen zposv cgemqrt ctgexc dppcon sgesdd_lwork stgsen_lwork zposvx cgeqp3 ctgsen dppsv sgesv stpmqrt zpotrf cgeqrf ctgsen_lwork dpptrf sgesvd stpqrt zpotri cgeqrf_lwork ctpmqrt dpptri sgesvd_lwork stpttf zpotrs cgeqrfp ctpqrt dpptrs sgesvx stpttr zppcon cgeqrfp_lwork ctpttf dpstf2 sgetc2 strsyl zppsv cgeqrt ctpttr dpstrf sgetrf strtri zpptrf cgerqf ctrsyl dpteqr sgetri strtrs zpptri cgesc2 ctrtri dptsv sgetri_lwork strttf zpptrs cgesdd ctrtrs dptsvx sgetrs strttp zpstf2 cgesdd_lwork ctrttf dpttrf sgges stzrzf zpstrf cgesv ctrttp dpttrs sggev stzrzf_lwork zpteqr cgesvd ctzrzf dsbev sgglse zgbsv zptsv cgesvd_lwork ctzrzf_lwork dsbevd sgglse_lwork zgbtrf zptsvx cgesvx cuncsd dsbevx sgtsv zgbtrs zpttrf cgetc2 cuncsd_lwork dsfrk sgtsvx zgebal zpttrs cgetrf cunghr dstebz sgttrf zgecon zrot cgetri cunghr_lwork dstein sgttrs zgeequ zsycon cgetri_lwork cungqr dstemr slamch zgeequb zsyconv cgetrs cungrq dstemr_lwork slange zgees zsyequb cgges cunmqr dsterf slarf zgeev zsysv cggev cunmrz dstev slarfg zgeev_lwork zsysv_lwork cgglse cunmrz_lwork dsycon slartg zgehrd zsysvx cgglse_lwork dgbsv dsyconv slasd4 zgehrd_lwork zsysvx_lwork cgtsv dgbtrf dsyequb slaswp zgels zsytf2 cgtsvx dgbtrs dsyev slauum zgels_lwork zsytrf cgttrf dgebal dsyev_lwork sorcsd zgelsd zsytrf_lwork cgttrs dgecon dsyevd sorcsd_lwork zgelsd_lwork ztbtrs chbevd dgeequ dsyevd_lwork sorghr zgelss ztfsm chbevx dgeequb dsyevr sorghr_lwork zgelss_lwork ztfttp checon dgees dsyevr_lwork sorgqr zgelsy ztfttr cheequb dgeev dsyevx sorgrq zgelsy_lwork ztgexc cheev dgeev_lwork dsyevx_lwork sormqr zgemqrt ztgsen cheev_lwork dgehrd dsygst sormrz zgeqp3 ztgsen_lwork cheevd dgehrd_lwork dsygv sormrz_lwork zgeqrf ztpmqrt cheevd_lwork dgejsv dsygv_lwork spbsv zgeqrf_lwork ztpqrt cheevr dgels dsygvd spbtrf zgeqrfp ztpttf cheevr_lwork dgels_lwork dsygvx spbtrs zgeqrfp_lwork ztpttr cheevx dgelsd dsygvx_lwork spftrf zgeqrt ztrsyl cheevx_lwork dgelsd_lwork dsysv spftri zgerqf ztrtri chegst dgelss dsysv_lwork spftrs zgesc2 ztrtrs chegv dgelss_lwork dsysvx spocon zgesdd ztrttf chegv_lwork dgelsy dsysvx_lwork sposv zgesdd_lwork ztrttp chegvd dgelsy_lwork dsytf2 sposvx zgesv ztzrzf chegvx dgemqrt dsytrd spotrf zgesvd ztzrzf_lwork chegvx_lwork dgeqp3 dsytrd_lwork spotri zgesvd_lwork zuncsd chesv dgeqrf dsytrf spotrs zgesvx zuncsd_lwork chesv_lwork dgeqrf_lwork dsytrf_lwork sppcon zgetc2 zunghr chesvx dgeqrfp dtbtrs sppsv zgetrf zunghr_lwork chesvx_lwork dgeqrfp_lwork dtfsm spptrf zgetri zungqr chetrd dgeqrt dtfttp spptri zgetri_lwork zungrq chetrd_lwork dgerqf dtfttr spptrs zgetrs zunmqr chetrf dgesc2 dtgexc spstf2 zgges zunmrz chetrf_lwork dgesdd dtgsen spstrf zggev zunmrz_lwork chfrk dgesdd_lwork dtgsen_lwork spteqr zgglse clange dgesv dtpmqrt sptsv zgglse_lwork clapack dgesvd dtpqrt sptsvx zgtsv clarf dgesvd_lwork dtpttf spttrf zgtsvx clarfg dgesvx dtpttr spttrs zgttrf
C++ API¶
Scipy wrapper API¶
Eigen¶
- https://eigen.tuxfamily.org/dox/
- Written in C++
- Previously used for linear algebra in Tensorflow (single CPU variant probably)
- Low and high-level (e.g. matrix decompositions) linear algebra functions
- Sparse linear algebra
- Can use BLAS/LaPACK as a backend https://eigen.tuxfamily.org/dox/TopicUsingBlasLapack.html
...name can be annoyingly google-proof however
Take-home message¶
Much (most) code you get for free, or even pay for, that implements mathematical algorithms may be of limited value.
Linear algebra, however, is "solved", in terms of numerical programming.
Your main concern is to find an implementation that fits for your application. E.g. parallel, or sparse, or optimized for a particular CPU.
Lab: Import and directly use packages in Python¶
- BLAS - compute a dot product of two random vectors
- LaPACK - solve a random linear system
II. Overview of leading Python tools¶
Leading Python tools:¶
Jupyter - "notebooks" for inline code + LaTex math + markup, etc.
NumPy - low-level array & matrix handling and algorithms
SciPy - higher level numerical algorithms (still fairly basic)
Matplotlib - matlab-style plotting & image display
SkLearn - (Scikit-learn) Machine Learning library
NumPy (short version)¶
- Numerical algorithm toolbox
- Similar to Matlab but Many key differences, such as zero-based indexing (like C) instead of one-based (like math texts).
- Uses BLAS/LAPACK typically
- NumPy Arrays - special data structure which allows direct and efficient linear algebra manipulations (basically a list with added functionality).
Harris, Charles R., et al. "Array programming with NumPy." nature 585.7825 (2020): 357-362.
Harris, Charles R., et al. "Array programming with NumPy." nature 585.7825 (2020): 357-362.
In [3]:
import numpy as np
np.show_config()
Build Dependencies:
blas:
detection method: pkgconfig
found: true
include directory: C:/Users/keith/anaconda3/envs/PREP_071325/Library/include
lib directory: C:/Users/keith/anaconda3/envs/PREP_071325/Library/lib
name: mkl-sdl
openblas configuration: unknown
pc file directory: C:\b\abs_958d_utj4g\croot\numpy_and_numpy_base_1750883811830\_h_env\Library\lib\pkgconfig
version: '2023.1'
lapack:
detection method: pkgconfig
found: true
include directory: C:/Users/keith/anaconda3/envs/PREP_071325/Library/include
lib directory: C:/Users/keith/anaconda3/envs/PREP_071325/Library/lib
name: mkl-sdl
openblas configuration: unknown
pc file directory: C:\b\abs_958d_utj4g\croot\numpy_and_numpy_base_1750883811830\_h_env\Library\lib\pkgconfig
version: '2023.1'
Compilers:
c:
commands: cl.exe
linker: link
name: msvc
version: 19.29.30159
c++:
commands: cl.exe
linker: link
name: msvc
version: 19.29.30159
cython:
commands: cython
linker: cython
name: cython
version: 3.0.11
Machine Information:
build:
cpu: x86_64
endian: little
family: x86_64
system: windows
host:
cpu: x86_64
endian: little
family: x86_64
system: windows
Python Information:
path: C:\b\abs_958d_utj4g\croot\numpy_and_numpy_base_1750883811830\_h_env\python.exe
version: '3.13'
SIMD Extensions:
baseline:
- SSE
- SSE2
- SSE3
found:
- SSSE3
- SSE41
- POPCNT
- SSE42
- AVX
- F16C
In [3]:
import numpy as np
x = np.array([2,0,1,8])
print(x)
[2 0 1 8]
Fast Numerical Mathematics¶
In [1]:
l = range(1234)
%timeit [i ** 2 for i in l]
389 µs ± 8.15 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
In [4]:
import numpy as np
a = np.arange(1234)
%timeit a ** 2
1.41 µs ± 7.46 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
NumPy for Matlab Users¶
In [4]:
help(np.array)
Help on built-in function array in module numpy.core.multiarray:
array(...)
array(object, dtype=None, copy=True, order='K', subok=False, ndmin=0)
Create an array.
Parameters
----------
object : array_like
An array, any object exposing the array interface, an object whose
__array__ method returns an array, or any (nested) sequence.
dtype : data-type, optional
The desired data-type for the array. If not given, then the type will
be determined as the minimum type required to hold the objects in the
sequence. This argument can only be used to 'upcast' the array. For
downcasting, use the .astype(t) method.
copy : bool, optional
If true (default), then the object is copied. Otherwise, a copy will
only be made if __array__ returns a copy, if obj is a nested sequence,
or if a copy is needed to satisfy any of the other requirements
(`dtype`, `order`, etc.).
order : {'K', 'A', 'C', 'F'}, optional
Specify the memory layout of the array. If object is not an array, the
newly created array will be in C order (row major) unless 'F' is
specified, in which case it will be in Fortran order (column major).
If object is an array the following holds.
===== ========= ===================================================
order no copy copy=True
===== ========= ===================================================
'K' unchanged F & C order preserved, otherwise most similar order
'A' unchanged F order if input is F and not C, otherwise C order
'C' C order C order
'F' F order F order
===== ========= ===================================================
When ``copy=False`` and a copy is made for other reasons, the result is
the same as if ``copy=True``, with some exceptions for `A`, see the
Notes section. The default order is 'K'.
subok : bool, optional
If True, then sub-classes will be passed-through, otherwise
the returned array will be forced to be a base-class array (default).
ndmin : int, optional
Specifies the minimum number of dimensions that the resulting
array should have. Ones will be pre-pended to the shape as
needed to meet this requirement.
Returns
-------
out : ndarray
An array object satisfying the specified requirements.
See Also
--------
empty, empty_like, zeros, zeros_like, ones, ones_like, full, full_like
Notes
-----
When order is 'A' and `object` is an array in neither 'C' nor 'F' order,
and a copy is forced by a change in dtype, then the order of the result is
not necessarily 'C' as expected. This is likely a bug.
Examples
--------
>>> np.array([1, 2, 3])
array([1, 2, 3])
Upcasting:
>>> np.array([1, 2, 3.0])
array([ 1., 2., 3.])
More than one dimension:
>>> np.array([[1, 2], [3, 4]])
array([[1, 2],
[3, 4]])
Minimum dimensions 2:
>>> np.array([1, 2, 3], ndmin=2)
array([[1, 2, 3]])
Type provided:
>>> np.array([1, 2, 3], dtype=complex)
array([ 1.+0.j, 2.+0.j, 3.+0.j])
Data-type consisting of more than one element:
>>> x = np.array([(1,2),(3,4)],dtype=[('a','<i4'),('b','<i4')])
>>> x['a']
array([1, 3])
Creating an array from sub-classes:
>>> np.array(np.mat('1 2; 3 4'))
array([[1, 2],
[3, 4]])
>>> np.array(np.mat('1 2; 3 4'), subok=True)
matrix([[1, 2],
[3, 4]])
Will cover more later, but first a quick intro to some other packages.
SciPy¶
Implements higher-level scientific algorithms using NumPy. Examples:
- Integration (scipy.integrate)
- Optimization (scipy.optimize)
- Interpolation (scipy.interpolate)
- Signal Processing (scipy.signal)
- Linear Algebra (scipy.linalg)
- Statistics (scipy.stats)
- File IO (scipy.io)
Also uses BLAS/LAPACK
In [30]:
import scipy
In [34]:
printcols(dir(scipy),2)
LowLevelCallable ndimage __numpy_version__ odr __version__ optimize cluster show_config fft signal fftpack sparse integrate spatial interpolate special io stats linalg test misc
In [35]:
printcols(dir(scipy.stats))
ConstantInputWarning bayes_mvs gausshyper matrix_normal rvs_ratio_uniforms DegenerateDataWarning bernoulli gaussian_kde maxwell scoreatpercentile FitError beta genexpon median_abs_deviation sem NearConstantInputWarning betabinom genextreme median_test semicircular NumericalInverseHermite betaprime gengamma mielke shapiro __all__ biasedurn genhalflogistic mode siegelslopes __builtins__ binned_statistic genhyperbolic moment sigmaclip __cached__ binned_statistic_2d geninvgauss monte_carlo_test skellam __doc__ binned_statistic_dd genlogistic mood skew __file__ binom gennorm morestats skewcauchy __loader__ binom_test genpareto moyal skewnorm __name__ binomtest geom mstats skewtest __package__ boltzmann gibrat mstats_basic somersd __path__ bootstrap gilbrat mstats_extras spearmanr __spec__ boschloo_exact gmean multinomial special_ortho_group _axis_nan_policy boxcox gompertz multiscale_graphcorr statlib _biasedurn boxcox_llf gstd multivariate_hypergeom stats _binned_statistic boxcox_normmax gumbel_l multivariate_normal studentized_range _binomtest boxcox_normplot gumbel_r multivariate_t t _boost bradford gzscore mvn test _common brunnermunzel halfcauchy mvsdist theilslopes _constants burr halfgennorm nakagami tiecorrect _continuous_distns burr12 halflogistic nbinom tmax _crosstab cauchy halfnorm ncf tmean _discrete_distns chi hmean nchypergeom_fisher tmin _distn_infrastructure chi2 hypergeom nchypergeom_wallenius trapezoid _distr_params chi2_contingency hypsecant nct trapz _entropy chisquare invgamma ncx2 triang _fit circmean invgauss nhypergeom trim1 _hypotests circstd invweibull norm trim_mean _hypotests_pythran circvar invwishart normaltest trimboth _kde combine_pvalues iqr norminvgauss truncexpon _ksstats contingency jarque_bera obrientransform truncnorm _levy_stable cosine johnsonsb ortho_group truncweibull_min _mannwhitneyu cramervonmises johnsonsu page_trend_test tsem _morestats cramervonmises_2samp kappa3 pareto tstd _mstats_basic crystalball kappa4 pearson3 ttest_1samp _mstats_extras cumfreq kde pearsonr ttest_ind _multivariate describe kendalltau percentileofscore ttest_ind_from_stats _mvn dgamma kruskal permutation_test ttest_rel _page_trend_test differential_entropy ks_1samp planck tukey_hsd _qmc dirichlet ks_2samp pmean tukeylambda _qmc_cy distributions ksone pointbiserialr tvar _relative_risk dlaplace kstat poisson uniform _resampling dweibull kstatvar power_divergence unitary_group _rvs_sampling energy_distance kstest powerlaw variation _sobol entropy kstwo powerlognorm vonmises _statlib epps_singleton_2samp kstwobign powernorm vonmises_line _stats erlang kurtosis ppcc_max wald _stats_mstats_common expon kurtosistest ppcc_plot wasserstein_distance _stats_py exponnorm laplace probplot weibull_max _tukeylambda_stats exponpow laplace_asymmetric qmc weibull_min _unuran exponweib levene randint weightedtau _variation f levy random_correlation wilcoxon _warnings_errors f_oneway levy_l rankdata wishart alexandergovern fatiguelife levy_stable ranksums wrapcauchy alpha find_repeats linregress rayleigh yeojohnson anderson fisher_exact loggamma rdist yeojohnson_llf anderson_ksamp fisk logistic recipinvgauss yeojohnson_normmax anglit fit loglaplace reciprocal yeojohnson_normplot ansari fligner lognorm relfreq yulesimon arcsine foldcauchy logser rice zipf argus foldnorm loguniform rv_continuous zipfian barnard_exact friedmanchisquare lomax rv_discrete zmap bartlett gamma mannwhitneyu rv_histogram zscore
Matplotlib¶
Tutorial from: https://github.com/amueller/scipy-2017-sklearn/blob/master/notebooks/02.Scientific_Computing_Tools_in_Python.ipynb
Another important part of machine learning is the visualization of data. The most common
tool for this in Python is matplotlib. It is an extremely flexible package, and
we will go over some basics here.
Jupyter has built-in "magic functions", the "matoplotlib inline" mode, which will draw the plots directly inside the notebook. Should be on by default.
In [5]:
%matplotlib inline
import matplotlib.pyplot as plt
# Plotting a line
x = np.linspace(0, 10, 100)
plt.plot(x, np.sin(x));
In [9]:
# Scatter-plot points
x = np.random.normal(size=500)
y = np.random.normal(size=500)
plt.scatter(x, y);
In [18]:
# Showing images using imshow
# - note that origin is at the top-left by default!
#x = np.linspace(1, 12, 100)
#y = x[:, np.newaxis]
#im = y * np.sin(x) * np.cos(y)
im = np.array([[1, 2, 3],[4,5,6],[6,7,8]])
import matplotlib.pyplot as plt
plt.imshow(im);
plt.colorbar();
plt.xlabel('x')
plt.ylabel('y')
plt.show();
In [19]:
# Contour plots
# - note that origin here is at the bottom-left by default!
plt.contour(im);
In [20]:
# 3D plotting
from mpl_toolkits.mplot3d import Axes3D
ax = plt.axes(projection='3d')
xgrid, ygrid = np.meshgrid(x, y.ravel())
ax.plot_surface(xgrid, ygrid, im, cmap=plt.cm.viridis, cstride=2, rstride=2, linewidth=0);
There are many more plot types available. See matplotlib gallery.
Test these examples: copy the Source Code
link, and put it in a notebook using the %load magic.
For example:
In [21]:
# %load http://matplotlib.org/mpl_examples/pylab_examples/ellipse_collection.py
SkLearn¶
Many Machine Learning functions.
Couple dozen core developers + hundreds of other contributors.
2011 tutorial has over 10,000 citations.
"Scikit-learn: Machine Learning in Python", Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, Jake Vanderplas, Alexandre Passos, David Cournapeau, Matthieu Brucher, Matthieu Perrot, Édouard Duchesnay; 12(Oct):2825−2830, 2011
Considered leading Machine Learning toolbox. Was getting discarded as field switched from single CPU to multicore, but appears to be making comeback with upgrades for scaling sizes and parallelism
In [41]:
# conda install anaconda::scikit-learn
import sklearn as sk
printcols(dir(sk),2)
__SKLEARN_SETUP__ base __all__ clone __builtins__ config_context __cached__ exceptions __check_build externals __doc__ get_config __file__ logger __loader__ logging __name__ os __package__ random __path__ set_config __spec__ setup_module __version__ show_versions _config sys _distributor_init utils
In [42]:
printcols(dir(sk.base),2)
BaseEstimator _check_feature_names_in BiclusterMixin _check_y ClassifierMixin _get_feature_names ClusterMixin _is_pairwise DensityMixin _num_features MetaEstimatorMixin _pprint MultiOutputMixin _safe_tags OutlierMixin check_X_y RegressorMixin check_array TransformerMixin clone _DEFAULT_TAGS copy _IS_32BIT defaultdict _OneToOneFeatureMixin estimator_html_repr _UnstableArchMixin get_config __builtins__ inspect __cached__ is_classifier __doc__ is_outlier_detector __file__ is_regressor __loader__ np __name__ platform __package__ re __spec__ warnings __version__
Machine Learning Framework using Sklearn¶
Given training data $(\mathbf x_{(i)},y_i)$ for $i=1,...,m$ --> lists of samples $\verb|X|$ and labels $\verb|y|$
Choose a model $f(\cdot)$ where we want to make $f(\mathbf x_{(i)})\approx y_i$ (for all $i$) --> choose sklearn estimator to use
Define a loss function $L(f(\mathbf x), y)$ to minimize by changing $f(\cdot)$ ...by adjusting the weights --> default choises for estimators, sometimes multiple options
The Sklearn API¶
sklearn has an Object Oriented interface
Most models/transforms/objects in sklearn are Estimator objects
In [2]:
class Estimator(object):
def fit(self, X, y=None):
"""Fit model to data X (and y)"""
self.some_attribute = self.some_fitting_method(X, y)
return self
def predict(self, X_test):
"""Make prediction based on passed features"""
pred = self.make_prediction(X_test)
return pred
model = Estimator()
The Estimator class defines a fit() method as well as a predict() method. For an instance of an Estimator stored in a variable model:
model.fit: fits the model with the passed in training data. For supervised models, it also accepts a second argumentythat corresponds to the labels (model.fit(X, y). For unsupervised models, there are no labels so you only need to pass in the feature matrix (model.fit(X))Since the interface is very OO, the instance itself stores the results of the
fitinternally. And as such you must alwaysfit()before youpredict()on the same object.model.predict: predicts new labels for any new datapoints passed in (model.predict(X_test)) and returns an array equal in length to the number of rows of what is passed in containing the predicted labels.
Types of subclass of estimator:¶
- Supervised
- Unsupervised
- Feature Processing
Supervised¶
Supervised estimators in addition to the above methods typically also have:
model.predict_proba: For classifiers that have a notion of probability (or some measure of confidence in a prediction) this method returns those "probabilities". The label with the highest probability is what is returned by themodel.predict()` mehod from above.model.score: For both classification and regression models, this method returns some measure of validation of the model (which is configurable). For example, in regression the default is typically R^2 and classification it is accuracy.
Unsupervised - Transformer interface¶
Some estimators in the library implement this.
Unsupervised in this case refers to any method that does not need labels, including unsupervised classifiers, preprocessing (like tf-idf), dimensionality reduction, etc.
The transformer interface usually defines two additional methods:
model.transform: Given an unsupervised model, transform the input into a new basis (or feature space). This accepts on argument (usually a feature matrix) and returns a matrix of the input transformed. Note: You need tofit()the model before you transform it.model.fit_transform: For some models you may not need tofit()andtransform()separately. In these cases it is more convenient to do both at the same time.
In [1]:
X = [[0], [1], [2], [3]]
y = [0, 0, 1, 1]
from sklearn.neighbors import KNeighborsClassifier
neigh = KNeighborsClassifier(n_neighbors=3)
neigh.fit(X, y)
print(neigh.predict([[1.1]]))
print(neigh.predict_proba([[0.9]]))
[0] [[0.66666667 0.33333333]]
Recap¶
- Mathematical software is implemented hierarchically
- BLAS = Lowest level(s), implemented on particular hardware, e.g. for x86 processors
- Higher level functions (solvers and decompositions) calls BLAS to implement internal computations
- Structured and Sparse matrix libraries also call BLAS typically (convert large sparse problems to small dense problems)
- Modern alternatives also can utilize BLAS/LAPACK backend and also support hardware-specific replacements (e.g. CUDA)
IV. NumPy Basics¶
NumPy Arrays for Vectors¶
In [23]:
x = np.array([2,0,1,8])
x
Out[23]:
array([2, 0, 1, 8])
NumPy arrays are zero-based like Python lists
In [5]:
print(x[0],x[3])
2 8
Matlab-ish slicing¶
In [6]:
A = np.array([[1,2,3],[4,5,6],[7,8,9]])
print(A)
[[1 2 3] [4 5 6] [7 8 9]]
In [7]:
print(A[:,0])
[1 4 7]
In [8]:
print(A[1,:])
[4 5 6]
In [10]:
print(A[1,0:2])
[4 5]
In [11]:
print(A[:,0:1])
[[1] [4] [7]]
Array programming with numpy¶
In [41]:
v = np.array([1,2,3])
w = np.array([3,4,5])
print(v)
print(5*v)
print(5*v+3*w)
[1 2 3] [ 5 10 15] [14 22 30]
$\verb|numpy.linalg|$¶
- Use BLAS and LAPACK to provide efficient low level implementations of standard linear algebra algorithms.
- Those libraries may be provided by NumPy itself using C versions of a subset of their reference implementations
- When possible, highly optimized libraries that take advantage of specialized processor functionality are preferred.
- Examples of such libraries are OpenBLAS, MKL (TM), and ATLAS.
https://docs.scipy.org/doc/numpy/reference/routines.linalg.html
Matrix and vector products¶
- dot(a, b[, out]) Dot product of two arrays.
- linalg.multi_dot(arrays) Compute the dot product of two or more arrays in a single function call, while automatically selecting the fastest evaluation order.
- vdot(a, b) Return the dot product of two vectors.
- inner(a, b) Inner product of two arrays.
- outer(a, b[, out]) Compute the outer product of two vectors.
- matmul(x1, x2, /[, out, casting, order, …]) Matrix product of two arrays.
- tensordot(a, b[, axes]) Compute tensor dot product along specified axes.
- einsum(subscripts, *operands[, out, dtype, …]) Evaluates the Einstein summation convention on the operands.
- einsum_path(subscripts, *operands[, optimize]) Evaluates the lowest cost contraction order for an einsum expression by considering the creation of intermediate arrays.
- linalg.matrix_power(a, n) Raise a square matrix to the (integer) power n.
- kron(a, b) Kronecker product of two arrays.
Decompositions¶
- linalg.cholesky(a) Cholesky decomposition.
- linalg.qr(a[, mode]) Compute the qr factorization of a matrix.
- linalg.svd(a[, full_matrices, compute_uv, …]) Singular Value Decomposition.
Matrix eigenvalues¶
- linalg.eig(a) Compute the eigenvalues and right eigenvectors of a square array.
- linalg.eigh(a[, UPLO]) Return the eigenvalues and eigenvectors of a complex Hermitian (conjugate symmetric) or a real symmetric matrix.
- linalg.eigvals(a) Compute the eigenvalues of a general matrix.
- linalg.eigvalsh(a[, UPLO]) Compute the eigenvalues of a complex Hermitian or real symmetric matrix.
Norms and other numbers¶
- linalg.norm(x[, ord, axis, keepdims]) Matrix or vector norm.
- linalg.cond(x[, p]) Compute the condition number of a matrix.
- linalg.det(a) Compute the determinant of an array.
- linalg.matrix_rank(M[, tol, hermitian]) Return matrix rank of array using SVD method
- linalg.slogdet(a) Compute the sign and (natural) logarithm of the determinant of an array.
- trace(a[, offset, axis1, axis2, dtype, out]) Return the sum along diagonals of the array.
Solving equations and inverting matrices¶
- linalg.solve(a, b) Solve a linear matrix equation, or system of linear scalar equations.
- linalg.tensorsolve(a, b[, axes]) Solve the tensor equation a x = b for x.
- linalg.lstsq(a, b[, rcond]) Return the least-squares solution to a linear matrix equation.
- linalg.inv(a) Compute the (multiplicative) inverse of a matrix.
- linalg.pinv(a[, rcond, hermitian]) Compute the (Moore-Penrose) pseudo-inverse of a matrix.
- linalg.tensorinv(a[, ind]) Compute the ‘inverse’ of an N-dimensional array.
"Vectorization"¶
Converting loops into linear algebra functions where loop is performed at lower level
E.g. instead of looping over a list to compute squares of elements, make into array and "square" array
In [7]:
v = [1,2,3,4,5]
v2 = []
for v_k in v:
v2.append(v_k**2)
v2
Out[7]:
[1, 4, 9, 16, 25]
In [8]:
np.array(v)**2
Out[8]:
array([ 1, 4, 9, 16, 25], dtype=int32)
In [14]:
v = range(0,1000)
%timeit v2 = list([x**2 for x in v])
266 µs ± 7.21 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
In [4]:
import numpy as np
v = np.arange(0,1000)
%timeit v2 = v**2
1.33 µs ± 18.4 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
Dot Product using NumPy¶
In [13]:
w = np.array([1, 2])
v = np.array([3, 4])
np.dot(w,v)
Out[13]:
11
In [14]:
w = np.array([0, 1])
v = np.array([1, 0])
np.dot(w,v)
Out[14]:
0
In [15]:
w = np.array([1, 2])
v = np.array([-2, 1])
np.dot(w,v)
Out[15]:
0
Matrix-vector multiplication using NumPy - CAUTION¶
NumPy implements a broadcast version of multiplication when using the "*" operator.
It guesses what you mean rather by seeking ways the shapes match, rather than giving error.
In [17]:
w = np.array([[2,-6],[-1, 4]])
v = np.array([1,0])
print(w)
print(w*v)
[[ 2 -6] [-1 4]] [[ 2 0] [-1 0]]
In [18]:
v+w
Out[18]:
array([[ 3, -6],
[ 0, 4]])
"Broadcasting"¶
Multiply matrices or vectrors in an element-wise manner by repeating "singular" dimensions.
Matlab and Numpy automatically do this now.
Also (and previously) in matlab, e.g.: "brcm(@plus, A, x). numpy.broadcast
For true "Matrix-vector" multiplication, use matmul or "@"¶
a special case of matrix-matrix multiplication
In [19]:
w = np.array([[ 2, -6],
[-1, 4]])
In [20]:
v = np.array([[+1],
[-1]])
In [21]:
np.matmul(w,v)
Out[21]:
array([[ 8],
[-5]])
In [56]:
w@v
Out[56]:
array([ 8, -5])
In [22]:
np.dot(w,v)
Out[22]:
array([[ 8],
[-5]])
Matrix-matrix multiplication¶
In [25]:
A = np.array([[1,2],[3, 4]])
B = np.array([[1,0],[1, 1]])
print(A*B)
[[1 0] [3 4]]
What is this doing?¶
In [6]:
np.matmul(A,B)
Out[6]:
array([[3, 2],
[7, 4]])
In [5]:
np.dot(A,B)
Out[5]:
array([[3, 2],
[7, 4]])
In [4]:
A@B
Out[4]:
array([[3, 2],
[7, 4]])
In [58]:
u = [1,2]
v = [1,3]
np.dot(np.array(u),np.array(v))
Out[58]:
7
Reshaping vectors into matrices and whatever else you want¶
Note also flatten() function to convert to "flat" array
In [23]:
A = np.array([[1,2,3],[4,5,6]])
print(A)
[[1 2 3] [4 5 6]]
In [27]:
A.flatten(), A.flatten().shape
Out[27]:
(array([1, 2, 3, 4, 5, 6]), (6,))
In [83]:
x = np.array([1,2,3,4])
x = x.reshape(4,1)
print(x.reshape(4,1))
[[1] [2] [3] [4]]
In [84]:
print(x.reshape(2,2))
[[1 2] [3 4]]
In [85]:
print(x)
x_T = x.transpose()
print (x_T)
print (x_T.shape)
print(x.T)
[[1] [2] [3] [4]] [[1 2 3 4]] (1, 4) [[1 2 3 4]]
In [48]:
x = np.array([1,2,3,4])
print(x.shape)
print(x)
print(x.T)
(4,) [1 2 3 4] [1 2 3 4]
In [47]:
x = np.array([[1,2,3,4]])
print(x.shape)
print(x)
print(x.T)
(1, 4) [[1 2 3 4]] [[1] [2] [3] [4]]
Famous matrices in NumPy¶
In [20]:
np.zeros((3,4))
Out[20]:
array([[0., 0., 0., 0.],
[0., 0., 0., 0.],
[0., 0., 0., 0.]])
In [21]:
np.ones((3,4))
Out[21]:
array([[1., 1., 1., 1.],
[1., 1., 1., 1.],
[1., 1., 1., 1.]])
In [16]:
np.eye(3)
Out[16]:
array([[1., 0., 0.],
[0., 1., 0.],
[0., 0., 1.]])
In [17]:
D = np.diag([1,2,3])
print(D)
[[1 0 0] [0 2 0] [0 0 3]]
In [18]:
np.diag(D)
Out[18]:
array([1, 2, 3])
In [22]:
# rand = uniform random numbers in [0,1]
np.random.rand(3,4) # note two inputs, not tuple
Out[22]:
array([[0.99147344, 0.01075122, 0.81974321, 0.10942128],
[0.94041317, 0.50408574, 0.71896232, 0.12113072],
[0.77483078, 0.32627551, 0.44489225, 0.13208849]])
In [23]:
# randn ~ N(0,1)
np.random.randn(3,4)
Out[23]:
array([[ 1.09642959, 1.40715723, -0.30932909, 0.66846735],
[ 0.22536484, -0.08283707, -0.60713286, -0.16836462],
[-1.20191307, 1.53390883, -0.48587343, -0.84379159]])
Lab: Numpy practice¶
Implement the "linear algebra with Python" Exercises from the bootcamp using numpy:
https://www.keithdillon.com/classes/UMMC/Bootcamp/__python_exercises_01.html