Программа трендов L1 не работает, CVXOPT, CVXPY - PullRequest
Новая
       24

Программа трендов L1 не работает, CVXOPT, CVXPY

0 голосов
/ 01 мая 2018

Я взял следующий код отсюда L1-фильтрация тренда

Теперь у меня есть Python 2.7, и мой код выглядит следующим образом: 

import cvxopt as cvxopt
import scipy as scipy
import scipy.sparse
import cvxpy as cvx
import numpy as np 
import matplotlib.pyplot as plt
import random

Я генерирую случайный сигнал здесь: 

amplitude = 10
t = 100
random.seed()
tau = random.uniform(3, 4)
X = np.arange(t)
noise = np.random.normal(0,1,100)
y = amplitude * np.exp(-X/tau)+noise

Следующий код, который я получил по ссылке. Этот код предназначен для оценки сигнала с использованием регуляризации l1-нормы, которая не работает для меня: 

n = y.size
# Form second difference matrix.
e = np.mat(np.ones((1, n)))
# D = scipy.sparse.spdiags(np.vstack((e, -2*e, e)), range(3), n-2, n)
D = scipy.sparse.spdiags(np.vstack((-e,e)), range(2), n-1, n)
# Convert D to cvxopt sparse format, due to bug in scipy which prevents
# overloading neccessary for CVXPY. Use COOrdinate format as intermediate.
D_coo = D.tocoo()
D = cvxopt.spmatrix(D_coo.data, D_coo.row.tolist(), D_coo.col.tolist())

# Set regularization parameter.
vlambda = 50
# Solve l1 trend filtering problem.
x = cvx.Variable(n)

obj = cvx.Minimize(0.5 * cvx.sum_squares(y - x) + vlambda*cvx.norm(D*x,1))
prob = cvx.Problem(obj)
# ECOS and SCS solvers fail to converge before
# the iteration limit. Use CVXOPT instead.
prob.solve(solver=cvx.CVXOPT,verbose=True)

print('Solver status: ', prob.status)
# Check for error.
if prob.status != cvx.OPTIMAL:
    raise Exception("Solver did not converge!")


# Plot estimated trend with original signal.
plt.figure(figsize=(6, 6))
plt.plot(np.arange(1,n+1), y, 'k:', linewidth=1.0)
plt.plot(np.arange(1,n+1), np.array(x.value), 'b-', linewidth=2.0)
plt.xlabel('date')
plt.ylabel('log price')
plt.show()

Я получаю следующую ошибку: 

    Traceback (most recent call last):
      File "MultiClipVideoEditing.py", line 81, in <module>
        obj = cvx.Minimize(0.5 * cvx.sum_squares(y - x) + vlambda*cvx.norm(D*x,1))
      File "C:\Users\klbm9\Anaconda3\envs\OptAndCv\lib\site-packages\cvxpy\expressions\expression.py", line 40, in cast_op
other = self.cast_to_const(other)
      File "C:\Users\klbm9\Anaconda3\envs\OptAndCv\lib\site-packages\cvxpy\expressions\expression.py", line 242, in cast_to_const
return expr if isinstance(expr, Expression) else cvxtypes.constant()(expr)
      File "C:\Users\klbm9\Anaconda3\envs\OptAndCv\lib\site-packages\cvxpy\expressions\constants\constant.py", line 39, in __init__
self._value = intf.DEFAULT_INTF.const_to_matrix(value)
      File "C:\Users\klbm9\Anaconda3\envs\OptAndCv\lib\site-packages\cvxpy\interface\base_matrix_interface.py", line 45, in new_converter
if not convert_scalars and 
    cvxpy.interface.matrix_utilities.is_scalar(value):
      File "C:\Users\klbm9\Anaconda3\envs\OptAndCv\lib\site-packages\cvxpy\interface\matrix_utilities.py", line 147, in is_scalar
return size(constant) == (1, 1)
      File "C:\Users\klbm9\Anaconda3\envs\OptAndCv\lib\site-packages\cvxpy\interface\matrix_utilities.py", line 135, in size
raise TypeError("%s is not a valid type for a Constant value." % type(constant))
    TypeError: <type 'cvxopt.base.spmatrix'> is not a valid type for a Constant value.

Я считаю, что при умножении D*x есть проблема, но я не знаю, как это исправить.

Ответы [ 2 ]

0 голосов
/ 11 мая 2018

Следующий код работает без проблем с CVXPY 1.0 и OSQP. Нет необходимости конвертировать какую-либо матрицу в плотный формат (это может быть очень-очень медленно!). Обратите внимание, что CVXPY 1.0 все еще проходит тестирование и в настоящее время работает только с Python 2.7. 

import scipy.sparse
import cvxpy as cvx
import numpy as np
import matplotlib.pyplot as plt
import random

# Generate random signal
amplitude = 10
t = 100
random.seed()
tau = random.uniform(3, 4)
X = np.arange(t)
noise = np.random.normal(0, 1, 100)
y = amplitude * np.exp(-X/tau)+noise
n = y.size

# Form second difference matrix.
e = np.mat(np.ones((1, n)))
D = scipy.sparse.spdiags(np.vstack((e, -2*e, e)), range(3), n-2, n)

# Set regularization parameter.
vlambda = 50

# Solve l1 trend filtering problem.
x = cvx.Variable(n)
obj = cvx.Minimize(0.5 * cvx.sum_squares(y - x) + vlambda*cvx.norm(D * x, 1))
prob = cvx.Problem(obj)

# Solve with OSQP (default solver)
prob.solve(verbose=True)

print('Solver status: ', prob.status)
# Check for error.
if prob.status != cvx.OPTIMAL:
    raise Exception("Solver did not converge!")


# Plot estimated trend with original signal.
plt.figure(figsize=(6, 6))
plt.plot(np.arange(1, n+1), y, 'k:', linewidth=1.0)
plt.plot(np.arange(1, n+1), np.array(x.value), 'b-', linewidth=2.0)
plt.xlabel('date')
plt.ylabel('log price')
plt.show()

См. график вывода

0 голосов
/ 01 мая 2018

Я нашел обходной путь, вместо того, чтобы использовать матрицу из cvxopt, т.е. 

    D = scipy.sparse.spdiags(np.vstack((-e,e)), range(2), n-1, n)
    D_coo = D.tocoo()
    D = cvxopt.spmatrix(D_coo.data, D_coo.row.tolist(), D_coo.col.tolist())

Я напрямую преобразовываю разреженную матрицу в плотную матрицу следующим образом: 

    from scipy.sparse import csr_matrix

    D = scipy.sparse.spdiags(np.vstack((-e,e)), range(2), n-1, n)
    D = csr_matrix(D)

Я использовал эту ссылку для преобразования скудной разреженной матрицы в плотную Это дает действительный D*x и решает мою проблему.

...