Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

simplify e_pri and e_dual computation using squared_norm properties a… #41

Open
wants to merge 3 commits into
base: master
Choose a base branch
from
Open

simplify e_pri and e_dual computation using squared_norm properties a… #41

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
9ccc218
simplify e_pri and e_dual computation using squared_norm properties a…
fdtomasi Jul 6, 2023
fae775b
add validation test
fdtomasi Jul 6, 2023
aec1adb
add math details on convergence criterion
fdtomasi Jul 11, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
115 changes: 46 additions & 69 deletions regain/covariance/kernel_latent_time_graphical_lasso_.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -153,32 +153,31 @@ def kernel_latent_time_graphical_lasso(
for m in range(1, n_times):
Z_L = Z_0.copy()[:-m]
Z_R = Z_0.copy()[m:]
Z_M[m] = (Z_L, Z_R)

W_L = np.zeros_like(Z_L)
W_R = np.zeros_like(Z_R)
W_M[m] = (W_L, W_R)
# Dims: 2 T-m D D
Z_M[m] = np.array((Z_L, Z_R))
W_M[m] = np.zeros_like(Z_M[m])

Y_L = np.zeros_like(Z_L)
Y_R = np.zeros_like(Z_R)
Y_M[m] = (Y_L, Y_R)
Y_M[m] = np.zeros_like(Z_M[m])
U_M[m] = np.zeros_like(W_M[m])

U_L = np.zeros_like(W_L)
U_R = np.zeros_like(W_R)
U_M[m] = (U_L, U_R)

Z_L_old = np.zeros_like(Z_L)
Z_R_old = np.zeros_like(Z_R)
Z_M_old[m] = (Z_L_old, Z_R_old)

W_L_old = np.zeros_like(W_L)
W_R_old = np.zeros_like(W_R)
W_M_old[m] = (W_L_old, W_R_old)
Z_M_old[m] = np.zeros_like(Z_M[m])
W_M_old[m] = np.zeros_like(W_M[m])

if n_samples is None:
n_samples = np.ones(n_times)

checks = []
# Convergence criterion for KLTGL.
# The number of auxiliary variables, Z_0, Z_M and U_M are:
# d^2 T + 2 d^2 T(T-1) / 2 + 2 d^2 T(T-1) / 2
# so:
# c = tol * [d^2 T + 4 d^2 T(T-1) / 2]^(1/2)
# = tol * d [T + 2 T(T-1)]^(1/2)
# = tol * d [T + 2T^2 - 2T]^(1/2)
# = tol * d [2T^2 - T]^(1/2)
# = tol * d [T (2T - 1)]^(1/2)
c = n_features * np.sqrt(n_times * (2 * n_times - 1)) * tol
for iteration_ in range(max_iter):
# update R
A = Z_0 - W_0 - X_0
Expand Down Expand Up @@ -237,7 +236,7 @@ def kernel_latent_time_graphical_lasso(
rtol=rtol,
max_iter=max_iter,
)
Z_M[m] = (Z_L, Z_R)
Z_M[m] = np.array((Z_L, Z_R))

# update other residuals
Y_L += Z_0[:-m] - Z_L
Expand All @@ -262,7 +261,7 @@ def kernel_latent_time_graphical_lasso(
rtol=rtol,
max_iter=max_iter,
)
W_M[m] = (W_L, W_R)
W_M[m] = np.array((W_L, W_R))

# update other residuals
U_L += W_0[:-m] - W_L
Expand All @@ -283,10 +282,7 @@ def kernel_latent_time_graphical_lasso(
snorm = rho * np.sqrt(
squared_norm(R - R_old)
+ sum(
squared_norm(Z_M[m][0] - Z_M_old[m][0])
+ squared_norm(Z_M[m][1] - Z_M_old[m][1])
+ squared_norm(W_M[m][0] - W_M_old[m][0])
+ squared_norm(W_M[m][1] - W_M_old[m][1])
squared_norm(Z_M[m] - Z_M_old[m]) + squared_norm(W_M[m] - W_M_old[m])
for m in range(1, n_times)
)
)
Expand All @@ -311,53 +307,39 @@ def kernel_latent_time_graphical_lasso(
else np.nan
)

check = Convergence(
obj=obj,
rnorm=rnorm,
snorm=snorm,
e_pri=n_features * np.sqrt(n_times * (2 * n_times - 1)) * tol
+ rtol
* max(
np.sqrt(
squared_norm(R)
+ sum(
squared_norm(Z_M[m][0])
+ squared_norm(Z_M[m][1])
+ squared_norm(W_M[m][0])
+ squared_norm(W_M[m][1])
for m in range(1, n_times)
)
),
np.sqrt(
squared_norm(Z_0 - W_0)
+ sum(
squared_norm(Z_0[:-m])
+ squared_norm(Z_0[m:])
+ squared_norm(W_0[:-m])
+ squared_norm(W_0[m:])
for m in range(1, n_times)
)
),
e_pri = c + rtol * max(
np.sqrt(
squared_norm(R)
+ sum(
squared_norm(Z_M[m]) + squared_norm(W_M[m])
for m in range(1, n_times)
)
),
e_dual=n_features * np.sqrt(n_times * (2 * n_times - 1)) * tol
+ rtol
* rho
* np.sqrt(
squared_norm(X_0)
np.sqrt(
squared_norm(Z_0 - W_0)
+ sum(
squared_norm(Y_M[m][0])
+ squared_norm(Y_M[m][1])
+ squared_norm(U_M[m][0])
+ squared_norm(U_M[m][1])
squared_norm(Z_0[:-m])
+ squared_norm(Z_0[m:])
+ squared_norm(W_0[:-m])
+ squared_norm(W_0[m:])
for m in range(1, n_times)
)
),
)
e_dual = c + rtol * rho * np.sqrt(
squared_norm(X_0)
+ sum(
squared_norm(Y_M[m]) + squared_norm(U_M[m]) for m in range(1, n_times)
)
)
check = Convergence(
obj=obj, rnorm=rnorm, snorm=snorm, e_pri=e_pri, e_dual=e_dual
)

R_old = R.copy()
for m in range(1, n_times):
Z_M_old[m] = (Z_M[m][0].copy(), Z_M[m][1].copy())
W_M_old[m] = (W_M[m][0].copy(), W_M[m][1].copy())
Z_M_old[m] = Z_M[m].copy()
W_M_old[m] = W_M[m].copy()

if verbose:
print(check)
Expand All @@ -372,13 +354,8 @@ def kernel_latent_time_graphical_lasso(
# scaled dual variables should be also rescaled
X_0 *= rho / rho_new
for m in range(1, n_times):
Y_L, Y_R = Y_M[m]
Y_L *= rho / rho_new
Y_R *= rho / rho_new

U_L, U_R = U_M[m]
U_L *= rho / rho_new
U_R *= rho / rho_new
Y_M[m] *= rho / rho_new
U_M[m] *= rho / rho_new
rho = rho_new
else:
warnings.warn("Objective did not converge.")
Expand Down
78 changes: 32 additions & 46 deletions regain/covariance/kernel_time_graphical_lasso_.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -149,15 +149,9 @@ def kernel_time_graphical_lasso(
# all possible markovians jumps
Z_L = Z_0.copy()[:-m]
Z_R = Z_0.copy()[m:]
Z_M[m] = (Z_L, Z_R)

U_L = np.zeros_like(Z_L)
U_R = np.zeros_like(Z_R)
U_M[m] = (U_L, U_R)

Z_L_old = np.zeros_like(Z_L)
Z_R_old = np.zeros_like(Z_R)
Z_M_old[m] = (Z_L_old, Z_R_old)
Z_M[m] = np.array((Z_L, Z_R))
U_M[m] = np.zeros_like(Z_M[m])
Z_M_old[m] = np.zeros_like(Z_M[m])

if n_samples is None:
n_samples = np.ones(n_times)
Expand All @@ -167,6 +161,17 @@ def kernel_time_graphical_lasso(
obj=objective(n_samples, emp_cov, Z_0, Z_0, Z_M, alpha, kernel, psi)
)
]

# Convergence criterion for KLTGL.
# The number of auxiliary variables, Z_0, Z_M are:
# d^2 T + 2 d^2 T(T-1) / 2
# so:
# c = tol * [d^2 T + 2 d^2 T(T-1) / 2]^(1/2)
# = tol * d [T + T(T-1)]^(1/2)
# = tol * d [T + T^2 - T]^(1/2)
# = tol * d [T^2]^(1/2)
# = tol * d T
c = n_features * n_times * tol
for iteration_ in range(max_iter):
# update K
A = Z_0 - U_0
Expand Down Expand Up @@ -215,7 +220,7 @@ def kernel_time_graphical_lasso(
rtol=rtol,
max_iter=max_iter,
)
Z_M[m] = (Z_L, Z_R)
Z_M[m] = np.array((Z_L, Z_R))

# update other residuals
U_L += K[:-m] - Z_L
Expand All @@ -232,11 +237,7 @@ def kernel_time_graphical_lasso(

snorm = rho * np.sqrt(
squared_norm(Z_0 - Z_0_old)
+ sum(
squared_norm(Z_M[m][0] - Z_M_old[m][0])
+ squared_norm(Z_M[m][1] - Z_M_old[m][1])
for m in range(1, n_times)
)
+ sum(squared_norm(Z_M[m] - Z_M_old[m]) for m in range(1, n_times))
)

obj = (
Expand All @@ -245,42 +246,29 @@ def kernel_time_graphical_lasso(
else np.nan
)

check = Convergence(
obj=obj,
rnorm=rnorm,
snorm=snorm,
e_pri=n_features * n_times * tol
+ rtol
* max(
np.sqrt(
squared_norm(Z_0)
+ sum(
squared_norm(Z_M[m][0]) + squared_norm(Z_M[m][1])
for m in range(1, n_times)
)
),
np.sqrt(
squared_norm(K)
+ sum(
squared_norm(K[:-m]) + squared_norm(K[m:])
for m in range(1, n_times)
)
),
# We use: squared_norm(Z_M[m]) == squared_norm(Z_M[m][0]) + squared_norm(Z_M[m][1])
e_pri = c + rtol * max(
np.sqrt(
squared_norm(Z_0) + sum(squared_norm(Z_M[m]) for m in range(1, n_times))
),
e_dual=n_features * n_times * tol
+ rtol
* rho
* np.sqrt(
squared_norm(U_0)
np.sqrt(
squared_norm(K)
+ sum(
squared_norm(U_M[m][0]) + squared_norm(U_M[m][1])
squared_norm(K[:-m]) + squared_norm(K[m:])
for m in range(1, n_times)
)
),
)
# We use: squared_norm(U_M[m]) == squared_norm(U_M[m][0]) + squared_norm(U_M[m][1])
e_dual = c + rtol * rho * np.sqrt(
squared_norm(U_0) + sum(squared_norm(U_M[m]) for m in range(1, n_times))
)
check = Convergence(
obj=obj, rnorm=rnorm, snorm=snorm, e_pri=e_pri, e_dual=e_dual
)
Z_0_old = Z_0.copy()
for m in range(1, n_times):
Z_M_old[m] = (Z_M[m][0].copy(), Z_M[m][1].copy())
Z_M_old[m] = Z_M[m].copy()

if verbose:
print(check)
Expand All @@ -299,9 +287,7 @@ def kernel_time_graphical_lasso(
# scaled dual variables should be also rescaled
U_0 *= rho / rho_new
for m in range(1, n_times):
U_L, U_R = U_M[m]
U_L *= rho / rho_new
U_R *= rho / rho_new
U_M[m] *= rho / rho_new
rho = rho_new
else:
warnings.warn("Objective did not converge.")
Expand Down
67 changes: 35 additions & 32 deletions regain/covariance/latent_time_graphical_lasso_.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -167,6 +167,16 @@ def latent_time_graphical_lasso(
n_samples = np.ones(emp_cov.shape[0])

checks = []
# Convergence criterion for LTGL.
# The number of auxiliary variables, Z_0, Z_1, Z_2, W_1, W_2 are:
# d^2 T + 4 d^2 (T-1)
# so
# c = tol * [d^2 T + 4 d^2 (T-1)]^(1/2)
# c = tol * d [T + 4 (T-1)]^(1/2)
# c = tol * d [T + 4T - 4]^(1/2)
# c = tol * d [5T - 4]^(1/2)
c = np.sqrt(Z_0.size + 4 * Z_1.size) * tol

for iteration_ in range(max_iter):
# update R
A = Z_0 - W_0 - X_0
Expand Down Expand Up @@ -277,41 +287,34 @@ def latent_time_graphical_lasso(
else np.nan
)

check = Convergence(
obj=obj,
rnorm=rnorm,
snorm=snorm,
e_pri=np.sqrt(R.size + 4 * Z_1.size) * tol
+ rtol
* max(
np.sqrt(
squared_norm(R)
+ squared_norm(Z_1)
+ squared_norm(Z_2)
+ squared_norm(W_1)
+ squared_norm(W_2)
),
np.sqrt(
squared_norm(Z_0 - W_0)
+ squared_norm(Z_0[:-1])
+ squared_norm(Z_0[1:])
+ squared_norm(W_0[:-1])
+ squared_norm(W_0[1:])
),
e_pri = c + rtol * max(
np.sqrt(
squared_norm(R)
+ squared_norm(Z_1)
+ squared_norm(Z_2)
+ squared_norm(W_1)
+ squared_norm(W_2)
),
e_dual=np.sqrt(R.size + 4 * Z_1.size) * tol
+ rtol
* rho
* (
np.sqrt(
squared_norm(X_0)
+ squared_norm(X_1)
+ squared_norm(X_2)
+ squared_norm(U_1)
+ squared_norm(U_2)
)
np.sqrt(
squared_norm(Z_0 - W_0)
+ squared_norm(Z_0[:-1])
+ squared_norm(Z_0[1:])
+ squared_norm(W_0[:-1])
+ squared_norm(W_0[1:])
),
)
e_dual = c + rtol * rho * (
np.sqrt(
squared_norm(X_0)
+ squared_norm(X_1)
+ squared_norm(X_2)
+ squared_norm(U_1)
+ squared_norm(U_2)
)
)
check = Convergence(
obj=obj, rnorm=rnorm, snorm=snorm, e_pri=e_pri, e_dual=e_dual
)

R_old = R.copy()
Z_1_old = Z_1.copy()
Expand Down
Loading