Revision

**241bf7ad2b806f6677a5e23534247f35f3a70f10**authored by rballester on**19 February 2023, 19:35:27 UTC**, committed by rballester on**19 February 2023, 19:35:27 UTC****1 parent**be80cb2

maxvol.py

```
# Copyright (c) 2015
# Alexander Mikhalev
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
# This code was shamelessly taken from https://bitbucket.org/muxas/maxvolpy
# in order to make the dependency optional. It is till recommended to install
# the original package as it has faster cython based implementations.
import numpy as np
from scipy.linalg import get_lapack_funcs, get_blas_funcs
def py_rect_maxvol(A, tol=1., maxK=None, min_add_K=None, minK=None,
start_maxvol_iters=10, identity_submatrix=True, top_k_index=-1):
"""
Python implementation of rectangular 2-volume maximization.
See Also
--------
rect_maxvol
"""
# tol2 - square of parameter tol
tol2 = tol**2
# N - number of rows, r - number of columns of matrix A
N, r = A.shape
# some work on parameters
if N <= r:
return np.arange(N, dtype=np.int32), np.eye(N, dtype=A.dtype)
if maxK is None or maxK > N:
maxK = N
if maxK < r:
maxK = r
if minK is None or minK < r:
minK = r
if minK > N:
minK = N
if min_add_K is not None:
minK = max(minK, r + min_add_K)
if minK > maxK:
minK = maxK
#raise ValueError('minK value cannot be greater than maxK value')
if top_k_index == -1 or top_k_index > N:
top_k_index = N
if top_k_index < r:
top_k_index = r
# choose initial submatrix and coefficients according to maxvol
# algorithm
index = np.zeros(N, dtype=np.int32)
chosen = np.ones(top_k_index)
tmp_index, C = py_maxvol(A, 1.05, start_maxvol_iters, top_k_index)
index[:r] = tmp_index
chosen[tmp_index] = 0
C = np.asfortranarray(C)
# compute square 2-norms of each row in coefficients matrix C
row_norm_sqr = np.array([chosen[i]*np.linalg.norm(C[i], 2)**2 for
i in range(top_k_index)])
# find maximum value in row_norm_sqr
i = np.argmax(row_norm_sqr)
K = r
# set cgeru or zgeru for complex numbers and dger or sger
# for float numbers
try:
ger = get_blas_funcs('geru', [C])
except:
ger = get_blas_funcs('ger', [C])
# augment maxvol submatrix with each iteration
while (row_norm_sqr[i] > tol2 and K < maxK) or K < minK:
# add i to index and recompute C and square norms of each row
# by SVM-formula
index[K] = i
chosen[i] = 0
c = C[i].copy()
v = C.dot(c.conj())
l = 1.0/(1+v[i])
ger(-l,v,c,a=C,overwrite_a=1)
C = np.hstack([C, l*v.reshape(-1,1)])
row_norm_sqr -= (l*v[:top_k_index]*v[:top_k_index].conj()).real
row_norm_sqr *= chosen
# find maximum value in row_norm_sqr
i = row_norm_sqr.argmax()
K += 1
# parameter identity_submatrix is True, set submatrix,
# corresponding to maxvol rows, equal to identity matrix
if identity_submatrix:
C[index[:K]] = np.eye(K, dtype=C.dtype)
return index[:K].copy(), C
def py_maxvol(A, tol=1.05, max_iters=100, top_k_index=-1):
"""
Python implementation of 1-volume maximization.
See Also
--------
maxvol
"""
# some work on parameters
if tol < 1:
tol = 1.0
N, r = A.shape
if N <= r:
return np.arange(N, dtype=np.int32), np.eye(N, dtype=A.dtype)
if top_k_index == -1 or top_k_index > N:
top_k_index = N
if top_k_index < r:
top_k_index = r
# set auxiliary matrices and get corresponding *GETRF function
# from lapack
B = np.copy(A[:top_k_index], order='F')
C = np.copy(A.T, order='F')
H, ipiv, info = get_lapack_funcs('getrf', [B])(B, overwrite_a=1)
# compute pivots from ipiv (result of *GETRF)
index = np.arange(N, dtype=np.int32)
for i in range(r):
tmp = index[i]
index[i] = index[ipiv[i]]
index[ipiv[i]] = tmp
# solve A = CH, H is in LU format
B = H[:r]
# It will be much faster to use *TRSM instead of *TRTRS
trtrs = get_lapack_funcs('trtrs', [B])
trtrs(B, C, trans=1, lower=0, unitdiag=0, overwrite_b=1)
trtrs(B, C, trans=1, lower=1, unitdiag=1, overwrite_b=1)
# C has shape (r, N) -- it is stored transposed
# find max value in C
i, j = divmod(abs(C[:,:top_k_index]).argmax(), top_k_index)
# set cgeru or zgeru for complex numbers and dger or sger for
# float numbers
try:
ger = get_blas_funcs('geru', [C])
except:
ger = get_blas_funcs('ger', [C])
# set number of iters to 0
iters = 0
# check if need to swap rows
while abs(C[i,j]) > tol and iters < max_iters:
# add j to index and recompute C by SVM-formula
index[i] = j
tmp_row = C[i].copy()
tmp_column = C[:,j].copy()
tmp_column[i] -= 1.
alpha = -1./C[i,j]
ger(alpha, tmp_column, tmp_row, a=C, overwrite_a=1)
iters += 1
i, j = divmod(abs(C[:,:top_k_index]).argmax(), top_k_index)
return index[:r].copy(), C.T
```

Computing file changes ...