A small Python package demonstrating how to use LAPACKE and CBLAS with NumPy arrays in C extension modules.
The project repository also includes an example of making LAPACKE calls from
normal C code in the lapacke_demo directory. See the corresponding
lapacke_demo/README.rst for details on building and running the example.
Important Note:
Although the code has been written to allow linking against Intel MKL, in practice this has not proven possible, as calls to Intel MKL routines result in segmentation faults. The linker line suggested by the Intel Link Line Advisor leads to fatal errors during runtime since some symbols cannot be found, while using the single dynamic library
libmkl_rt.soappears to [occasionally] not segfault whenMKL_INTERFACE_LAYERis set toGNU,ILP64orGNU,LP64. The dedicated Intel article gives further details on setting the Intel MKL interface and threading layer.
TBD. However, the C extension modules can be built with either OpenBLAS or
standard system CBLAS and LAPACKE implementations. But note that unless linked
against Intel MKL using ILP64 interface [1] or with OpenBLAS using 64-bit
int, i.e. built with INTERFACE64=1, no npypacke function or method
that accepts NumPy arrays should be passed arrays requiring 64-bit indexing as
an OverflowError will be raised. Do note that this package is for
demonstration purposes and so shouldn't be used for true "big data"
applications anyways.
| [1] | The ILP64 interface for Intel MKL uses 64-bit integers. See the Intel article on ILP64 vs. LP64 for more details. |
TBD. Local x86-64 builds on WSL Ubuntu 20.04 LTS linking against OpenBLAS
0.3.17 or against libblas 3.9.0 and liblapacke 3.9.0 installed using
apt have succeeded, as have builds on GitHub Actions runners for
manylinux1, manylinux2010, and MacOS linked against the latest OpenBLAS
(0.3.17), built using cibuildwheel. However, there have been issues with
building wheels for Windows with the latest OpenBLAS Windows binaries.
The npypacke package contains the regression and solvers
subpackages. The regression subpackage provides the LinearRegression
class, implemented like a scikit-learn estimator, which can be used to fit
a linear model with optional intercept by ordinary least squares, using either
QR or singular value decompositions. The solvers subpackage provides the
mnewton local minimizer, implemented such that it may be used as a frontend
for scipy.optimize.minimize by passing mnewton to the method
keyword argument. mnewton implements Newton's method with a Hessian
modification, where any Hessians that are not positive definite have a multiple
of the identity added to them to make the Newton direction using the modified
Hessian a descent direction. mnewton implements Algorithm 3.3 on page 51 of
Nocedal and Wright's Numerical Optimization and uses the returned lower
Cholesky factor of the modified Hessian when computing the descent direction.
The step size for the line search satisfies the Armijo condition and is chosen
using a backtracking line search.
Here are some examples that demonstrate how the LinearRegression class and
mnewton minimizer could be used. To run the first example, scikit-learn
must be installed. NumPy and SciPy must of course be installed.
Fit a linear model with intercept by least squares on the Boston house prices data using QR decomposition.
from sklearn.datasets import load_boston
from sklearn.model_selection import train_test_split
from npypacke.regression import LinearRegression
X, y = load_boston(return_X_y=True)
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=7)
lr = LinearRegression(solver="qr").fit(X_train, y_train)
print(lr.score(X_test, y_test))Minimize the multivariate Rosenbrock function using Newton's method with diagonal Hessian modification.
import numpy as np
from scipy.optimize import rosen, rosen_der, rosen_hess
from npypacke.solvers import mnewton
res = mnewton(rosen, np.zeros(5), jac=rosen_der, hess=rosen_hess)
print(res.x)