1

I took the following code from here L1 trend filtering

Now i have python 2.7, and my code is as following :

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

I am generating a random signal here :

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

The following code i got from the link. This code is for estimating the signal using l1-norm regularisation, which is not working for me:

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()

I get the following error :

    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.

I believe there is a problem while multiplying D*x, but i dont know how to fix that.

Community
  • 1
  • 1
klbm9999
  • 31
  • 1
  • 9

2 Answers2

1

The following code works with no problems with CVXPY 1.0 and OSQP. No need to convert any matrix to dense format (it can be very very slow!). Note that CVXPY 1.0 is still under testing and currently works only with 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()

See the output plot

bstellato
  • 160
  • 1
  • 6
  • Hey, sorry for the late reply, but your solution really worked with a little edit!! and its way more faster than converting to a dense matrix (45mins~ 2secs, i am working with different input). Your code didnt work for me (error: cvx doesnt have an attribute OSQP), according to [this](http://osqp.readthedocs.io/en/latest/parsers/cvxpy.html) link, cvxpy has OSQP as the default solver (ie change to `prob.solve(verbose=True)`). Your answer was great help, thanks, please edit i will accept the answer. Also, i made a mistake, i am working with python 2.7, i made that edit too. – klbm9999 May 29 '18 at 07:38
  • Note that now CVXPY 1.0 is out now. I suggest you to install it. It is also working with Python 3.x. CVXPY 1.0 has the attribute OSQP because it is the default solver. If you got an error it means you were using an older version of CVXPY. Disclamer: I am one of the OSQP developers :) – bstellato May 30 '18 at 14:01
  • Thanks for that, will definitely check out. That helped a lot. – klbm9999 May 30 '18 at 14:25
0

I have found a workaround, rather than using the matrix from cvxopt, ie

    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())

I directly conver the sparse matrix to a dense matrix in this manner :

    from scipy.sparse import csr_matrix

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

I used this reference to convert a scipy sparse matrix to a dense one This gives a valid D*x and solves my problem.

klbm9999
  • 31
  • 1
  • 9