///|
fn training_decision(
index : Int,
labels : Array[Double],
alphas : Array[Double],
matrix : Array[Array[Double]],
bias : Double,
) -> Double {
let mut total = bias
for support = 0; support < alphas.length(); support = support + 1 {
if alphas[support] != 0.0 {
total = total + alphas[support] * labels[support] * matrix[support][index]
}
}
total
}
///|
fn choose_second_index(
first : Int,
first_error : Double,
labels : Array[Double],
alphas : Array[Double],
matrix : Array[Array[Double]],
bias : Double,
) -> Int {
let mut best = if first == 0 { 1 } else { 0 }
let mut best_difference = -1.0
for candidate = 0; candidate < labels.length(); candidate = candidate + 1 {
if candidate != first {
let candidate_error = training_decision(
candidate, labels, alphas, matrix, bias,
) -
labels[candidate]
let difference = (first_error - candidate_error).abs()
if difference > best_difference {
best_difference = difference
best = candidate
}
}
}
best
}
///|
fn clipped(value : Double, lower : Double, upper : Double) -> Double {
if value < lower {
lower
} else if value > upper {
upper
} else {
value
}
}
///|
fn maximum_kkt_violation(
labels : Array[Double],
alphas : Array[Double],
bounds : Array[Double],
matrix : Array[Array[Double]],
bias : Double,
epsilon : Double,
) -> Double {
let mut maximum = 0.0
for index = 0; index < labels.length(); index = index + 1 {
let margin = labels[index] *
training_decision(index, labels, alphas, matrix, bias)
let violation = if alphas[index] <= epsilon {
(1.0 - margin).max(0.0)
} else if alphas[index] >= bounds[index] - epsilon {
(margin - 1.0).max(0.0)
} else {
(margin - 1.0).abs()
}
maximum = maximum.max(violation)
}
maximum
}
///|
fn dual_objective(
labels : Array[Double],
alphas : Array[Double],
matrix : Array[Array[Double]],
) -> Double {
let mut linear = 0.0
let mut quadratic = 0.0
for left = 0; left < alphas.length(); left = left + 1 {
linear = linear + alphas[left]
for right = 0; right < alphas.length(); right = right + 1 {
quadratic = quadratic +
alphas[left] *
alphas[right] *
labels[left] *
labels[right] *
matrix[left][right]
}
}
linear - 0.5 * quadratic
}
///|
fn train_binary_bounds(
data : Dataset,
config : BinaryConfig,
bounds : Array[Double],
) -> Result[BinaryModel, SvmError] {
let classes = data.classes()
if classes.length() != 2 {
return Err(InvalidBinaryClassCount(classes.length()))
}
let rows = data.features()
let original_labels = data.labels()
let labels = Array::make(data.row_count(), -1.0)
for index, label in original_labels {
if label == classes[1] {
labels[index] = 1.0
}
}
let matrix = match kernel_matrix(config.kernel(), data) {
Err(error) => return Err(error)
Ok(value) => value
}
let alphas = Array::make(data.row_count(), 0.0)
let bias_box = [0.0]
let mut iterations = 0
let mut unchanged_passes = 0
let mut pair_updates = 0
while iterations < config.max_iterations() &&
unchanged_passes < config.max_passes() {
let mut changed = 0
for first = 0; first < data.row_count(); first = first + 1 {
let first_error = training_decision(
first,
labels,
alphas,
matrix,
bias_box[0],
) -
labels[first]
let violates = (
labels[first] * first_error < -config.tolerance() &&
alphas[first] < bounds[first] - config.alpha_epsilon()
) ||
(
labels[first] * first_error > config.tolerance() &&
alphas[first] > config.alpha_epsilon()
)
if violates {
let second = choose_second_index(
first,
first_error,
labels,
alphas,
matrix,
bias_box[0],
)
let second_error = training_decision(
second,
labels,
alphas,
matrix,
bias_box[0],
) -
labels[second]
let old_first = alphas[first]
let old_second = alphas[second]
let same_class = labels[first] == labels[second]
let (lower, upper) = if !same_class {
(
(old_second - old_first).max(0.0),
bounds[second].min(bounds[first] + old_second - old_first),
)
} else {
(
(old_first + old_second - bounds[first]).max(0.0),
bounds[second].min(old_first + old_second),
)
}
if upper > lower {
let eta = 2.0 * matrix[first][second] -
matrix[first][first] -
matrix[second][second]
if eta < 0.0 {
let proposed_second = old_second -
labels[second] * (first_error - second_error) / eta
let new_second = clipped(proposed_second, lower, upper)
if (new_second - old_second).abs() >
config.alpha_epsilon() *
(new_second + old_second + config.alpha_epsilon()) {
let new_first = old_first +
labels[first] * labels[second] * (old_second - new_second)
alphas[first] = clipped(new_first, 0.0, bounds[first])
alphas[second] = new_second
let first_delta = alphas[first] - old_first
let second_delta = alphas[second] - old_second
let first_bias = bias_box[0] -
first_error -
labels[first] * first_delta * matrix[first][first] -
labels[second] * second_delta * matrix[first][second]
let second_bias = bias_box[0] -
second_error -
labels[first] * first_delta * matrix[first][second] -
labels[second] * second_delta * matrix[second][second]
bias_box[0] = if alphas[first] > config.alpha_epsilon() &&
alphas[first] < bounds[first] - config.alpha_epsilon() {
first_bias
} else if alphas[second] > config.alpha_epsilon() &&
alphas[second] < bounds[second] - config.alpha_epsilon() {
second_bias
} else {
(first_bias + second_bias) / 2.0
}
changed = changed + 1
pair_updates = pair_updates + 1
}
}
}
}
}
iterations = iterations + 1
if changed == 0 {
unchanged_passes = unchanged_passes + 1
} else {
unchanged_passes = 0
}
}
if !finite_double(bias_box[0]) {
return Err(NumericFailure("SMO bias update"))
}
let support_rows : Array[Array[Double]] = []
let support_labels : Array[Int] = []
let support_alphas : Array[Double] = []
let coefficients : Array[Double] = []
for index, alpha in alphas {
if alpha > config.alpha_epsilon() {
support_rows.push(rows[index].copy())
support_labels.push(original_labels[index])
support_alphas.push(alpha)
coefficients.push(alpha * labels[index])
}
}
let maximum_violation = maximum_kkt_violation(
labels,
alphas,
bounds,
matrix,
bias_box[0],
config.alpha_epsilon(),
)
let objective = dual_objective(labels, alphas, matrix)
if !finite_double(maximum_violation) || !finite_double(objective) {
return Err(NumericFailure("SMO convergence report"))
}
let converged = unchanged_passes >= config.max_passes() ||
maximum_violation <= config.tolerance()
Ok({
support_rows,
support_class_labels: support_labels,
support_alpha_values: support_alphas,
signed_coefficients: coefficients,
bias_value: bias_box[0],
input_columns: data.feature_count(),
model_kernel: config.kernel(),
lower_class: classes[0],
upper_class: classes[1],
convergence_data: {
completed_iterations: iterations,
completed_pair_updates: pair_updates,
retained_support_vectors: support_rows.length(),
largest_kkt_violation: maximum_violation,
final_dual_objective: objective,
reached_stopping_condition: converged,
},
})
}
///|
/// Fits a deterministic, unweighted binary C-SVC model.
pub fn train_binary(
data : Dataset,
config : BinaryConfig,
) -> Result[BinaryModel, SvmError] {
let bounds = Array::make(data.row_count(), config.c())
train_binary_bounds(data, config, bounds)
}