(основано на примерах к CVXPY, ссылка)
In [3]:
import numpy as np
import cvxpy as cvx
print(cvx.installed_solvers())
import matplotlib.pyplot as plt
In [4]:
USE_COLAB = False
if USE_COLAB:
!wget -q https://github.com/amkatrutsa/MIPT-Opt/blob/master/01-Intro/kandinskiy_comp_viii.jpg?raw=true -O kandinskiy_comp_viii.jpg
!wget -q https://github.com/amkatrutsa/MIPT-Opt/blob/master/01-Intro/lena_small.jpg?raw=true -O lena_small.jpg
!wget -q https://github.com/amkatrutsa/MIPT-Opt/blob/master/01-Intro/test_sound.wav?raw=true -O test_sound.wav
else:
plt.rc("text", usetex=True)
In [5]:
from PIL import Image
# np.random.seed(1)
# Load image
orig_img = Image.open("./kandinskiy_comp_viii.jpg")
# orig_img = Image.open("./lena_small.jpg")
# Convert to arrays
Uorig = np.array(orig_img)
rows, cols, colors = Uorig.shape
print(Uorig.shape)
# Known is 1 if the pixel is known,
# 0 if the pixel was corrupted
Known = np.zeros((rows, cols, colors))
for i in range(rows):
for j in range(cols):
if np.random.random() > 0.7:
for k in range(colors):
Known[i, j, k] = 1
# Create corrupted image
Ucorr = Known * Uorig
# Convert to image
corr_img = Image.fromarray(np.uint8(Ucorr))
In [6]:
fig, ax = plt.subplots(1, 2, figsize=(20, 20))
ax[0].imshow(orig_img);
ax[0].set_title("Original image", fontsize=24)
ax[0].axis('off')
ax[1].imshow(corr_img);
ax[1].set_title("Corrupted image", fontsize=24)
ax[1].axis('off');
In [4]:
# Define variables and constraints
variables = []
constraints = []
for i in range(colors):
# Slice of image tensor for the i-th color
U = cvx.Variable(shape=(rows, cols))
# Append to the list of variables
variables.append(U)
# Append constraints: known pixels in the given image
# has to be equal to the pixels in reconstructed image
constraints.append(cvx.multiply(Known[:, :, i], U) == \
cvx.multiply(Known[:, :, i], Ucorr[:, :, i]))
In [6]:
# Create minimization problem with created constraints and solve it
prob = cvx.Problem(cvx.Minimize(cvx.tv(*variables)), constraints)
prob.solve(verbose=True, solver="SCS", max_iters=500)
Out[6]:
In [7]:
# Load variables values into array
rec_arr = np.zeros((rows, cols, colors), dtype=np.uint8)
for i in range(colors):
rec_arr[:, :, i] = variables[i].value
In [8]:
fig, ax = plt.subplots(1, 2, figsize=(15, 10))
img_rec = Image.fromarray(rec_arr)
ax[0].imshow(img_rec);
ax[0].set_title("In-Painted image", fontsize=24)
ax[0].axis('off')
img_diff = Image.fromarray(np.abs(Uorig - rec_arr))
ax[1].imshow(img_diff);
ax[1].set_title("Difference image", fontsize=24)
ax[1].axis('off');
In [10]:
import scipy.io.wavfile as siowav
import IPython.display as ipd
# Load original track
sr, sound = siowav.read("./test_sound.wav")
print(sr, sound.shape, type(sound[0]))
start = 500000
fin = 1000000
sound = sound.astype("float32")
# Cut and normalize original track
sound = sound[start:fin] / np.max(sound[start:fin])
n = sound.shape[0]
print("Problem dimension = {}".format(n))
plt.plot(sound)
# plt.xscale("log")
Out[10]:
In [11]:
ipd.Audio(sound, rate=sr)
Out[11]:
In [12]:
# Create corrupted signal and normalize it
corrupt_sound = sound + 0.1 * np.random.randn(sound.shape[0])
corrupt_sound = corrupt_sound / np.max(np.abs(corrupt_sound))
plt.plot(corrupt_sound)
# plt.xscale("log")
Out[12]:
In [13]:
ipd.Audio(corrupt_sound, rate=sr)
Out[13]:
In [24]:
# Create variable for reconstructed signal
x_rec = cvx.Variable(corrupt_sound.shape[0])
# Create ac problem and solve it, \lambda = 1
prob = cvx.Problem(cvx.Minimize(cvx.norm(x_rec - corrupt_sound) + cvx.tv(x_rec)),
[x_rec >= -1, x_rec <= 1])
prob.solve(verbose=True, solver=cvx.SCS, max_iters=500)
Out[24]:
In [25]:
plt.plot(x_rec.value)
Out[25]:
In [26]:
ipd.Audio(x_rec.value, rate=sr)
Out[26]:
In [27]:
# The same as in the previous cell, but with \ell_2 regularization
x_rec = cvx.Variable(n)
prob = cvx.Problem(cvx.Minimize(cvx.norm(x_rec - corrupt_sound) +
cvx.square(cvx.norm(x_rec[1:n] - x_rec[0:n-1]))),
[x_rec >= -1, x_rec <= 1])
prob.solve(verbose=True, solver="SCS", max_iters=500)
Out[27]:
In [28]:
plt.plot(x_rec.value)
Out[28]:
In [29]:
ipd.Audio(x_rec.value, rate=sr)
Out[29]:
где $v_{t} = (v^x_t, v^y_t)$ и $p_t = (p^x_t, p^y_t)$
Приближённо имеем следующую систему
$$ v_{t+1} = \left(1 - \frac{\tau}{m}\eta \right)v_t + \frac{\tau}{m}u_t \quad p_{t+1} = p_t + \tau v_t $$Для записи в стандартной форме
$$ x_{t} = \begin{bmatrix} p^x_t \\ p^y_t \\ v^x_t \\ v^y_t \end{bmatrix} \quad A = \begin{bmatrix} 1 & 0 & \tau & 0 \\ 0 & 1 & 0 & \tau \\ 0 & 0 & 1 - \frac{\tau}{m}\eta & 0\\ 0 & 0 & 0 & 1 - \frac{\tau}{m}\eta \end{bmatrix} \quad B = \begin{bmatrix} 0 & 0 \\ 0 & 0 \\ \frac{\tau}{m} & 0 \\ 0 & \frac{\tau}{m} \end{bmatrix} $$
In [19]:
m = 1
eta = 0.1
tau = 0.01
T = 100
A = np.array([[1, 0, tau, 0],
[0, 1, 0, tau],
[0, 0, 1 - eta * tau / m, 0],
[0, 0, 0, 1 - eta * tau / m]])
B = np.array([[0, 0],
[0, 0],
[tau / m, 0],
[0, tau / m]])
In [20]:
# Initialize original and target states of the defined system
x0 = np.array([5, 5, 20, -5])
x_target = np.array([0, 7, 0, 0])
In [21]:
# Create variables to store evolution of the system and
# all control vectors
x = cvx.Variable((4, T + 1))
u = cvx.Variable((2, T))
states = []
for t in range(T):
# Create objective function and constraints in every time step
cost = cvx.sum_squares(u[:, t])
constr = [x[:, t+1] == A*x[:, t] + B*u[:, t]]
states.append(cvx.Problem(cvx.Minimize(cost), constr))
# Sum-up all problems: sum objectives and concatenate constraints
prob0 = sum(states)
# Append original state and target state
complete_constr = prob0.constraints + [x[:,0] == x0, x[:, T] == x_target]
prob = cvx.Problem(prob0.objective, complete_constr)
prob.solve(verbose=True)
Out[21]:
In [22]:
plt.plot(x.value[0, :], x.value[1, :])
plt.grid(True)
plt.xticks(fontsize=24)
plt.yticks(fontsize=24)
plt.xlabel(r"$x$", fontsize=24)
plt.ylabel(r"$y$", fontsize=24)
Out[22]:
In [23]:
plt.plot(u.value[0, :], label=r"$u_x$")
plt.plot(u.value[1, :], label=r"$u_y$")
plt.legend(fontsize=24)
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.xlabel("Time, $t$", fontsize=24)
Out[23]:
In [24]:
m1 = 1
eta1 = 0.1
m2 = 1
eta2 = 1
tau = 0.01
T = 100
A = np.array([[1, 0, tau, 0],
[0, 1, 0, tau],
[0, 0, 1 - eta1 * tau / m1, 0],
[0, 0, 0, 1 - eta1 * tau / m1]])
B = np.array([[0, 0],
[0, 0],
[tau / m1, 0],
[0, tau / m1]])
C = np.array([[1, 0, tau, 0],
[0, 1, 0, tau],
[0, 0, 1 - eta2 * tau / m2, 0],
[0, 0, 0, 1 - eta2 * tau / m2]])
D = np.array([[0, 0],
[0, 0],
[tau / m2, 0],
[0, tau / m2]])
In [25]:
x0 = np.array([0, 0, 1, 1])
z0 = np.array([10, 0, -1, 1])
In [26]:
x = cvx.Variable((4, T + 1))
z = cvx.Variable((4, T + 1))
u = cvx.Variable((2, T))
v = cvx.Variable((2, T))
states = []
for t in range(T):
cost = cvx.sum_squares(u[:, t]) + cvx.sum_squares(v[:, t])
constr = [x[:, t+1] == A*x[:, t] + B*u[:, t], z[:, t+1] == C*z[:, t] + D*v[:, t]]
states.append(cvx.Problem(cvx.Minimize(cost), constr))
prob0 = sum(states)
complete_constr = prob0.constraints + [x[:, T] == z[:, T], x[:, 0] == x0, z[:, 0] == z0]
prob = cvx.Problem(prob0.objective, complete_constr)
prob.solve(verbose=True)
Out[26]:
In [27]:
plt.plot(x.value[0, :], x.value[1, :], label="1")
plt.plot(z.value[0, :], z.value[1, :], label="2")
plt.legend(fontsize=24)
plt.xlabel(r"$x$", fontsize=24)
plt.ylabel(r"$y$", fontsize=24)
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
Out[27]:
Задача: дан набор точек на плоскости $x_i \in \mathbb{R}^2$, необходимо найти эллипс минимальной площади, такой что все точки $x_i$ лежат внутри него.
Приложения
Параметры:
эллипс можно задать аффинным преобразованием
$$ \{x \; | \; \| x \|^2_2 \leq 1\} \to \{ u \;| \; \|u\|^2_2 \leq 1, \; u = Ax + b\} $$Тогда площадь увеличивается в $\det (A^{-1})$ раз.
In [2]:
num_points = 50
n = 2
X = 2 * np.random.randn(n, num_points)
plt.figure(figsize=(10, 7))
plt.scatter(X[0, :], X[1, :], s=90)
Out[2]:
In [3]:
# Create variables: positive semi-definite matrix A and vector b
A = cvx.Variable((n,n), PSD=True)
b = cvx.Variable((n,))
prob = cvx.Problem(cvx.Minimize(-cvx.log_det(A)),
[cvx.norm(A * X[:, i] + b) <= 1 for i in range(num_points)])
In [4]:
prob.solve(verbose=True)
Out[4]:
In [6]:
plt.figure(figsize=(10, 7))
plt.scatter(X[0, :], X[1, :])
phi = np.linspace(0, 2 * np.pi, num=100)
xy = np.vstack((np.cos(phi) - b.value[0], np.sin(phi) - b.value[1]))
ellips = np.linalg.solve(A.value, xy)
plt.plot(ellips[0, :], ellips[1, :])
Out[6]:
In [61]:
import sklearn.datasets as skldata
import sklearn.model_selection as skms
np.random.seed(8)
n = 2
m = 2000
X, y = skldata.make_blobs(n_features=n, n_samples=m, centers=2, cluster_std=4)
y = 2 * y - 1
X_train, X_test, y_train, y_test = skms.train_test_split(X, y, test_size=0.5)
plt.scatter(X_train[:, 0], X_train[:, 1], c=y_train)
m_train = m // 2
m_test = m_train
In [62]:
beta = cvx.Variable(n)
v = cvx.Variable()
loss = cvx.sum(cvx.pos(1 - cvx.multiply(y_train, X_train*beta + v)))
reg = cvx.norm(beta, 2)
lambd = cvx.Parameter(nonneg=True, value=1)
prob = cvx.Problem(cvx.Minimize(loss / m_train + lambd * reg))
prob.solve(verbose=True)
Out[62]:
In [63]:
x = np.linspace(min(X_train[:, 0]), max(X_train[:, 0]))
plt.plot(x, (-v.value - beta.value[0] * x) / beta.value[1])
plt.scatter(X_train[:, 0], X_train[:, 1], c=y_train)
Out[63]:
In [42]:
import sklearn.datasets as skdata
import sklearn.preprocessing as skprep
X, y = skdata.load_breast_cancer(return_X_y=True)
X = skprep.StandardScaler().fit_transform(X)
y = 2 * y - 1
print("Data shape = {}".format(X.shape))
X_train, X_test, y_train, y_test = skms.train_test_split(X, y, test_size=0.3)
m_train, n = X_train.shape
m_test = X_test.shape[0]
# Redefine the problem
beta = cvx.Variable(n)
v = cvx.Variable()
loss = cvx.sum(cvx.pos(1 - cvx.multiply(y_train, X_train*beta + v)))
reg = cvx.norm(beta, 1)
lambd = cvx.Parameter(nonneg=True)
prob = cvx.Problem(cvx.Minimize(loss / m_train + lambd * reg))
In [52]:
TRIALS = 100
train_error = np.zeros(TRIALS)
test_error = np.zeros(TRIALS)
lambda_vals = np.logspace(-4, -1, TRIALS)
beta_vals = []
for i in range(TRIALS):
lambd.value = lambda_vals[i]
prob.solve(solver="ECOS")
train_error[i] = (y_train != np.sign(X_train.dot(beta.value) + v.value)).sum() * 1.0 / m_train
test_error[i] = (y_test != np.sign(X_test.dot(beta.value) + v.value)).sum() * 1.0 / m_test
beta_vals.append(beta.value)
In [53]:
plt.figure(figsize=(10, 7))
plt.plot(lambda_vals, train_error, label="Train error")
plt.plot(lambda_vals, test_error, label="Test error")
plt.xscale('log')
plt.legend(loc='best', fontsize=24)
plt.xlabel(r"$\lambda$", fontsize=24)
plt.xticks(fontsize=24)
_ = plt.yticks(fontsize=24)
In [54]:
plt.figure(figsize=(10, 7))
for i in range(n):
plt.plot(lambda_vals, [wi[i] for wi in beta_vals])
plt.xlabel(r"$\lambda$", fontsize=24)
plt.ylabel(r"$w_i$", fontsize=24)
plt.xscale("log")
plt.xticks(fontsize=24)
_ = plt.yticks(fontsize=24)