python解带L1正则的最小二乘

给定 $H\in {\mathbb{R}}^{d\times n}$，要解：$$\begin{gathered} \underset{Z}{\min }\parallel H-HZ{\parallel }_F^2+\alpha \parallel Z{\parallel }_1 \\ s.t.\text{diag}(Z)=0 \end{gathered}$$ 拆开每一列单独求解，变成：$$\underset{z_i}{\min }\parallel h_i-\hat{H}{z}_i{\parallel }_2^2+\alpha \parallel z_i{\parallel }_1$$ 其中 $\hat{H}$ 是把 H 第 i 列 $h_i$ 换成全 0。即带 L1 范数正则的最小二乘（L1-norm regularized least-squares）^[1,2]，可以用 cvxopt 包配合 [3] 实现好的 l1regls 求；或者用 [10]，是 pytorch 的方案。

Preliminary

安装见 [4]，conda、pip 都行。

mutual conversion

numpy.ndarray 和 cvxopt.matrix 之间的相互转换，参考 [5,6,7]。

import cvxopt
import numpy as np

H_np = np.arange(12).reshape(3, 4)
print(H_np)

# np.ndarray -> cvxopt.matrix
H_cvx = cvxopt.matrix(H_np)
print(H_cvx)

# cvxopt.matrix -> np.ndarray
print(np.array(H_cvx))

Code

cvxopt

• 数据量大时很慢，不支持权重 $\alpha$。
• 注意 H 要排成列向量。
• float 型 array 的 dtype 要转成 np.double 再转 cvxopt.matrix，见 [7]。
• [3] 要单独下载，且在开头加上 solvers.options['show_progress'] = False 让 solver 求解时不要 bb^[9]。

import struct
import cvxopt
import numpy as np
import l1regls  # downloaded from [3]


"""L1 regularized least square
min || H - HZ ||_2 + || Z ||_1
"""


# column vector
H = np.random.permutation(12).reshape(2, 6).astype(np.double)  # use `double`
d, n = H.shape
print("H:\n", H)

Z = np.zeros([n, n], dtype=np.double)
for i in range(n):
    # the vector to approximate
    h = H[:, i]

    # mask out the target to avoid trivial solution
    _mask = np.ones_like(H[0:1])
    _mask[0][i] = 0
    H_masked = H * _mask
    # print(H_masked)

    # np.ndarray -> cvxopt.matrix
    H_masked = cvxopt.matrix(H_masked)
    h = cvxopt.matrix(h)

    # solve
    z = l1regls.l1regls(H_masked, h)

    # cvxopt.matrix -> np.ndarray
    z = np.array(z)
    # z = np.vectorize(lambda x: struct.unpack('d', x))(_tmp)
    # z = np.vectorize(lambda x: int.from_bytes(x, 'big'))(np.array(z).T)

    Z[:, i] = z.flatten()


print("Z:\n", Z)
print("recon H:\n", H.dot(Z))

• 输出，可以看到，HZ 基本能重构出 H。

H:
[[ 6. 11.  4.  7.  8.  1.]
 [10.  0.  9.  3.  2.  5.]]
 
Z:
[[ 0.00000000e+00 -6.98943781e-06  8.92271576e-01  2.97727171e-01
   1.97727272e-01  1.59974460e-06]
 [ 1.38700201e-01  0.00000000e+00 -1.18924004e-01  4.69834622e-01
   6.15289256e-01 -1.03918158e-01]
 [ 1.10717194e+00 -3.29324447e-01  0.00000000e+00  6.28937736e-08
   4.00154515e-10  5.47136286e-01]
 [ 1.66159302e-05  1.62159367e-05 -4.83780134e-07  0.00000000e+00
   1.43251121e-09 -1.72704896e-08]
 [ 4.66304188e-06  1.52891479e+00 -1.34402595e-06  1.63529549e-07
   0.00000000e+00 -4.41351100e-08]
 [ 7.52874506e-06 -2.84235943e-05  3.14914688e-06  2.54879563e-08
   3.46453896e-10  0.00000000e+00]]

recon H:
[[ 5.95455111 10.91406371  4.04545442  6.95454546  7.95454545  1.04545454]
 [ 9.96464426  0.0937462   8.92272736  2.97727274  1.97727273  4.92424243]]

rfeinman/pytorch-lasso

• [10] 是 pytorch 实现的，支持 gpu、支持权重 $\alpha$。
• 如果用 pytorch 1.2，可能会报错说 ModuleNotFoundError: No module named 'torch._vmap_internals'，参考 [11]，把 lasso/__init__.py 中注掉 nonlinear 和 conv2d，这里只用到 linear。
• sparse_encode 第一个参数是 $h_i$，但支持传整个 H 进去，排成行向量；第二个参数是字典，本文即 $\hat{H}$，排成列向量，求得的 z 又是行向量，所以最后的重构公式是 $\hat{H}=ZH$。

import torch
from lasso.linear import sparse_encode#, dict_learning


n, d = 6, 2
# data, as ROW vector
H = torch.arange(n * d).view(n, d).float().cuda()
print("data:\n", H)

# Dictionary Learning
# dictionary, losses = dict_learning(data, n_components=50, alpha=0.5, algorithm='ista',
#   device='cuda', progbar=False)
# print(dictionary.size())

# dictionary, as COLUMN vector
D = H.T
# coefficient matrix, as ROW vector
Z = torch.zeros(n, n, dtype=torch.float).cuda()

for i in range(n):
    # the vector to approximate
    h = H[i:i+1]  # row vector

    # mask out the target to avoid trivial solution
    _mask = torch.ones_like(D[0:1])
    _mask[0][i] = 0
    D_masked = D * _mask
    # print(D_masked)

    z = sparse_encode(h, D_masked, alpha=0.2, algorithm='interior-point')
    # print(z.size())
    Z[i, :] = z

print("Z:\n", Z)

# H^ = ZH, instead of HZ
H_hat = Z.mm(H)
print("recon:\n", H_hat)

References

1. L1-norm regularized least-squares
2. L1-norm regularized least-squares on Python
3. l1regls.py
4. Installation instructions
5. Numpy and CVXOPT
6. python3 conversion between cvxopt.matrix and numpy.array
7. converting numpy vector to cvxopt
8. Fast l1-minimization algorithms and an application in robust face recognition: A review
9. How to silent cvxopt solver [Python]?
10. rfeinman/pytorch-lasso
11. ModuleNotFoundError: No module named ‘torch._vmap_internals’ with pytorch 1.2 #1
12. torch.lstsq