[Python知识库]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 求。

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

• 注意 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]]

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. torch.lstsq