Skip to content

Commit

Permalink
Performance improvement
Browse files Browse the repository at this point in the history
  • Loading branch information
@doccstat
doccstat committed Apr 23, 2024
1 parent 225565c commit 1918953
Show file tree
Hide file tree
Showing 6 changed files with 440 additions and 296 deletions.
321 changes: 321 additions & 0 deletions src/fastcpd_class_nll.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,327 @@

namespace fastcpd::classes {

colvec Fastcpd::get_gradient_arma(
const unsigned int segment_start,
const unsigned int segment_end,
const colvec& theta
) {
const mat data_segment = data.rows(segment_start, segment_end);
const unsigned int segment_length = segment_end - segment_start + 1;
mat reversed_data = reverse(data_segment, 0);
colvec reversed_theta = reverse(theta);
if (segment_length < max(order) + 1) {
return ones(theta.n_elem);
}
colvec variance_term = zeros(segment_length);
for (unsigned int i = max(order); i < segment_length; i++) {
variance_term(i) = data_segment(i, 0) - dot(
reversed_theta.rows(order(1) + 1, sum(order)),
data_segment.rows(i - order(0), i - 1)
) - dot(
reversed_theta.rows(1, order(1)),
variance_term.rows(i - order(1), i - 1)
);
}
colvec reversed_variance_term = reverse(variance_term);
mat phi_coefficient = zeros(segment_length, order(0)),
psi_coefficient = zeros(segment_length, order(1));
for (unsigned int i = max(order); i < segment_length; i++) {
phi_coefficient.row(i) = -reversed_data.rows(
segment_length - i, segment_length - i + order(0) - 1
).t() - reversed_theta.rows(1, order(1)).t() *
phi_coefficient.rows(i - order(1), i - 1);
}
for (unsigned int i = order(1); i < segment_length; i++) {
psi_coefficient.row(i) = -reversed_variance_term.rows(
segment_length - i, segment_length - i + order(1) - 1
).t() - reversed_theta.rows(1, order(1)).t() *
psi_coefficient.rows(i - order(1), i - 1);
}
colvec gradient = zeros(sum(order) + 1);
gradient.rows(0, order(0) - 1) =
phi_coefficient.row(segment_length - 1).t() *
variance_term(segment_length - 1) / theta(sum(order));
gradient.rows(order(0), sum(order) - 1) =
psi_coefficient.row(segment_length - 1).t() *
variance_term(segment_length - 1) / theta(sum(order));
gradient(sum(order)) = 1.0 / 2.0 / theta(sum(order)) -
std::pow(variance_term(segment_length - 1), 2) / 2.0 /
std::pow(theta(sum(order)), 2);
return gradient;
}

colvec Fastcpd::get_gradient_binomial(
const unsigned int segment_start,
const unsigned int segment_end,
const colvec& theta
) {
const mat data_segment = data.rows(segment_start, segment_end);
const unsigned int segment_length = segment_end - segment_start + 1;
rowvec new_data = data_segment.row(segment_length - 1);
rowvec x = new_data.tail(new_data.n_elem - 1);
double y = new_data(0);
return - (y - 1 / (1 + exp(-as_scalar(x * theta)))) * x.t();
}

colvec Fastcpd::get_gradient_lm(
const unsigned int segment_start,
const unsigned int segment_end,
const colvec& theta
) {
const mat data_segment = data.rows(segment_start, segment_end);
const unsigned int segment_length = segment_end - segment_start + 1;
rowvec new_data = data_segment.row(segment_length - 1);
rowvec x = new_data.tail(new_data.n_elem - 1);
double y = new_data(0);
return - (y - as_scalar(x * theta)) * x.t();
}

colvec Fastcpd::get_gradient_ma(
const unsigned int segment_start,
const unsigned int segment_end,
const colvec& theta
) {
const mat data_segment = data.rows(segment_start, segment_end);
const unsigned int segment_length = segment_end - segment_start + 1;
const unsigned int q = order(1);
mat reversed_data = reverse(data_segment, 0);
colvec reversed_theta = reverse(theta);
if (segment_length < q + 1) {
return ones(theta.n_elem);
}
colvec variance_term = zeros(segment_length);
for (unsigned int i = q; i < segment_length; i++) {
variance_term(i) = data_segment(i, 0) -
dot(reversed_theta.rows(1, q), variance_term.rows(i - q, i - 1));
}
colvec reversed_variance_term = reverse(variance_term);
mat psi_coefficient = zeros(segment_length, q);
for (unsigned int i = q; i < segment_length; i++) {
psi_coefficient.row(i) = -reversed_variance_term.rows(
segment_length - i, segment_length - i + q - 1
).t() - reversed_theta.rows(1, q).t() *
psi_coefficient.rows(i - q, i - 1);
}
colvec gradient = zeros(q + 1);
gradient.rows(0, q - 1) =
psi_coefficient.row(segment_length - 1).t() *
variance_term(segment_length - 1) / theta(q);
gradient(q) = 1.0 / 2.0 / theta(q) -
std::pow(variance_term(segment_length - 1), 2) / 2.0 /
std::pow(theta(q), 2);
return gradient;
}

colvec Fastcpd::get_gradient_poisson(
const unsigned int segment_start,
const unsigned int segment_end,
const colvec& theta
) {
const mat data_segment = data.rows(segment_start, segment_end);
const unsigned int segment_length = segment_end - segment_start + 1;
rowvec new_data = data_segment.row(segment_length - 1);
rowvec x = new_data.tail(new_data.n_elem - 1);
double y = new_data(0);
return - (y - exp(as_scalar(x * theta))) * x.t();
}

mat Fastcpd::get_hessian_arma(
const unsigned int segment_start,
const unsigned int segment_end,
const colvec& theta
) {
const mat data_segment = data.rows(segment_start, segment_end);
const unsigned int segment_length = segment_end - segment_start + 1;
// TODO(doccstat): Maybe we can store all these computations
mat reversed_data = reverse(data_segment, 0);
colvec reversed_theta = reverse(theta);
if (segment_length < max(order) + 1) {
return eye(theta.n_elem, theta.n_elem);
}
colvec variance_term = zeros(segment_length);
for (unsigned int i = max(order); i < segment_length; i++) {
variance_term(i) = data_segment(i, 0) - dot(
reversed_theta.rows(order(1) + 1, sum(order)),
data_segment.rows(i - order(0), i - 1)
) - dot(
reversed_theta.rows(1, order(1)),
variance_term.rows(i - order(1), i - 1)
);
}
colvec reversed_variance_term = reverse(variance_term);
mat phi_coefficient = zeros(segment_length, order(0)),
psi_coefficient = zeros(segment_length, order(1));
for (unsigned int i = max(order); i < segment_length; i++) {
phi_coefficient.row(i) = -reversed_data.rows(
segment_length - i, segment_length - i + order(0) - 1
).t() - reversed_theta.rows(1, order(1)).t() *
phi_coefficient.rows(i - order(1), i - 1);
}
for (unsigned int i = order(1); i < segment_length; i++) {
psi_coefficient.row(i) = -reversed_variance_term.rows(
segment_length - i, segment_length - i + order(1) - 1
).t() - reversed_theta.rows(1, order(1)).t() *
psi_coefficient.rows(i - order(1), i - 1);
}
mat reversed_coef_phi = reverse(phi_coefficient, 0),
reversed_coef_psi = reverse(psi_coefficient, 0);
cube phi_psi_coefficient = zeros(order(1), order(0), segment_length),
psi_psi_coefficient = zeros(order(1), order(1), segment_length);
for (unsigned int i = order(1); i < segment_length; i++) {
mat phi_psi_coefficient_part = zeros(order(1), order(0)),
psi_psi_coefficient_part = zeros(order(1), order(1));
for (unsigned int j = 1; j <= order(1); j++) {
phi_psi_coefficient_part +=
phi_psi_coefficient.slice(i - j) * theta(order(0) - 1 + j);
}
phi_psi_coefficient.slice(i) = -reversed_coef_phi.rows(
segment_length - i, segment_length - i + order(1) - 1
) - phi_psi_coefficient_part;
for (unsigned int j = 1; j <= order(1); j++) {
psi_psi_coefficient_part +=
psi_psi_coefficient.slice(i - j) * theta(order(0) - 1 + j);
}
psi_psi_coefficient.slice(i) = -reversed_coef_psi.rows(
segment_length - i, segment_length - i + order(1) - 1
) - reversed_coef_psi.rows(
segment_length - i, segment_length - i + order(1) - 1
).t() - psi_psi_coefficient_part;
}
mat hessian = zeros(sum(order) + 1, sum(order) + 1);
hessian.submat(0, 0, order(0) - 1, order(0) - 1) =
phi_coefficient.row(segment_length - 1).t() *
phi_coefficient.row(segment_length - 1) / theta(sum(order));
hessian.submat(0, order(0), order(0) - 1, sum(order) - 1) = (
phi_psi_coefficient.slice(segment_length - 1).t() *
variance_term(segment_length - 1) +
phi_coefficient.row(segment_length - 1).t() *
psi_coefficient.row(segment_length - 1)
) / theta(sum(order));
hessian.submat(order(0), 0, sum(order) - 1, order(0) - 1) =
hessian.submat(0, order(0), order(0) - 1, sum(order) - 1).t();
hessian.submat(0, sum(order), order(0) - 1, sum(order)) =
-phi_coefficient.row(segment_length - 1).t() *
variance_term(segment_length - 1) / theta(sum(order)) / theta(sum(order));
hessian.submat(sum(order), 0, sum(order), order(0) - 1) =
hessian.submat(0, sum(order), order(0) - 1, sum(order)).t();
hessian.submat(order(0), order(0), sum(order) - 1, sum(order) - 1) = (
psi_coefficient.row(segment_length - 1).t() *
psi_coefficient.row(segment_length - 1) +
psi_psi_coefficient.slice(segment_length - 1) *
variance_term(segment_length - 1)
) / theta(sum(order));
hessian.submat(order(0), sum(order), sum(order) - 1, sum(order)) =
-psi_coefficient.row(segment_length - 1).t() *
variance_term(segment_length - 1) / theta(sum(order)) / theta(sum(order));
hessian.submat(sum(order), order(0), sum(order), sum(order) - 1) =
hessian.submat(order(0), sum(order), sum(order) - 1, sum(order)).t();
hessian(sum(order), sum(order)) =
std::pow(variance_term(segment_length - 1), 2) /
std::pow(theta(sum(order)), 3) -
1.0 / 2.0 / std::pow(theta(sum(order)), 2);
return hessian;
}

mat Fastcpd::get_hessian_binomial(
const unsigned int segment_start,
const unsigned int segment_end,
const colvec& theta
) {
const mat data_segment = data.rows(segment_start, segment_end);
const unsigned int segment_length = segment_end - segment_start + 1;
rowvec new_data = data_segment.row(segment_length - 1);
rowvec x = new_data.tail(new_data.n_elem - 1);
double prob = 1 / (1 + exp(-dot(x, theta)));
return (x.t() * x) * ((1 - prob) * prob);
}

mat Fastcpd::get_hessian_lm(
const unsigned int segment_start,
const unsigned int segment_end,
const colvec& theta
) {
const mat data_segment = data.rows(segment_start, segment_end);
const unsigned int segment_length = segment_end - segment_start + 1;
rowvec new_data = data_segment.row(segment_length - 1);
rowvec x = new_data.tail(new_data.n_elem - 1);
return x.t() * x;
}

mat Fastcpd::get_hessian_ma(
const unsigned int segment_start,
const unsigned int segment_end,
const colvec& theta
) {
const mat data_segment = data.rows(segment_start, segment_end);
const unsigned int segment_length = segment_end - segment_start + 1;
const unsigned int q = order(1);
// TODO(doccstat): Maybe we can store all these computations
mat reversed_data = reverse(data_segment, 0);
colvec reversed_theta = reverse(theta);
if (segment_length < q + 1) {
return eye(theta.n_elem, theta.n_elem);
}
colvec variance_term = zeros(segment_length);
for (unsigned int i = q; i < segment_length; i++) {
variance_term(i) = data_segment(i, 0) - dot(
reversed_theta.rows(1, q), variance_term.rows(i - q, i - 1)
);
}
colvec reversed_variance_term = reverse(variance_term);
mat psi_coefficient = zeros(segment_length, q);
for (unsigned int i = q; i < segment_length; i++) {
psi_coefficient.row(i) = -reversed_variance_term.rows(
segment_length - i, segment_length - i + q - 1
).t() - reversed_theta.rows(1, q).t() *
psi_coefficient.rows(i - q, i - 1);
}
mat reversed_coef_psi = reverse(psi_coefficient, 0);
cube psi_psi_coefficient = zeros(q, q, segment_length);
for (unsigned int i = q; i < segment_length; i++) {
mat psi_psi_coefficient_part = zeros(q, q);
for (unsigned int j = 1; j <= q; j++) {
psi_psi_coefficient_part +=
psi_psi_coefficient.slice(i - j) * theta(j - 1);
}
psi_psi_coefficient.slice(i) = -reversed_coef_psi.rows(
segment_length - i, segment_length - i + q - 1
) - reversed_coef_psi.rows(
segment_length - i, segment_length - i + q - 1
).t() - psi_psi_coefficient_part;
}
mat hessian = zeros(q + 1, q + 1);
hessian.submat(0, 0, q - 1, q - 1) = (
psi_coefficient.row(segment_length - 1).t() *
psi_coefficient.row(segment_length - 1) +
psi_psi_coefficient.slice(segment_length - 1) *
variance_term(segment_length - 1)
) / theta(q);
hessian.submat(0, q, q - 1, q) =
-psi_coefficient.row(segment_length - 1).t() *
variance_term(segment_length - 1) / theta(q) / theta(q);
hessian.submat(q, 0, q, q - 1) = hessian.submat(0, q, q - 1, q).t();
hessian(q, q) =
std::pow(variance_term(segment_length - 1), 2) /
std::pow(theta(q), 3) -
1.0 / 2.0 / std::pow(theta(q), 2);
return hessian;
}

mat Fastcpd::get_hessian_poisson(
const unsigned int segment_start,
const unsigned int segment_end,
const colvec& theta
) {
const mat data_segment = data.rows(segment_start, segment_end);
const unsigned int segment_length = segment_end - segment_start + 1;
rowvec new_data = data_segment.row(segment_length - 1);
rowvec x = new_data.tail(new_data.n_elem - 1);
double prob = exp(as_scalar(x * theta));
// Prevent numerical issues if `prob` is too large.
return (x.t() * x) * std::min(as_scalar(prob), 1e10);
}

CostResult Fastcpd::get_nll_arma(
const unsigned int segment_start,
const unsigned int segment_end
Expand Down
Loading

0 comments on commit 1918953

Please sign in to comment.