diff --git a/include/equilibration.h b/include/equilibration.h index 9a595abe..9a831f76 100644 --- a/include/equilibration.h +++ b/include/equilibration.h @@ -41,6 +41,16 @@ void ruiz_equilibration(QOCOProblemData* data, QOCOScaling* scaling, QOCOInt ruiz_iters); +/** + * @brief Applies row normalization to scale data matrices. Sets E and F so + * that each row of A and each inequality row of G has unit infinity norm. + * Does not modify D or the cost scaling factor k. + * + * @param data Pointer to problem data. + * @param scaling Pointer to scaling struct. + */ +void row_normalization(QOCOProblemData* data, QOCOScaling* scaling); + /** * @brief Undo variable transformation induced by ruiz equilibration. * diff --git a/src/equilibration.c b/src/equilibration.c index b02d8bd5..72afb98f 100644 --- a/src/equilibration.c +++ b/src/equilibration.c @@ -179,6 +179,89 @@ void ruiz_equilibration(QOCOProblemData* data, QOCOScaling* scaling, sync_vector_to_device(scaling->Finvruiz); } +void row_normalization(QOCOProblemData* data, QOCOScaling* scaling) +{ + set_cpu_mode(1); + + // Initialize all scaling to identity. + for (QOCOInt i = 0; i < data->n; ++i) { + set_element_vectorf(scaling->Druiz, i, 1.0); + set_element_vectorf(scaling->Dinvruiz, i, 1.0); + } + for (QOCOInt i = 0; i < data->p; ++i) { + set_element_vectorf(scaling->Eruiz, i, 1.0); + set_element_vectorf(scaling->Einvruiz, i, 1.0); + } + for (QOCOInt i = 0; i < data->m; ++i) { + set_element_vectorf(scaling->Fruiz, i, 1.0); + set_element_vectorf(scaling->Finvruiz, i, 1.0); + } + scaling->k = 1.0; + scaling->kinv = 1.0; + + QOCOFloat* Eruiz_data = get_data_vectorf(scaling->Eruiz); + QOCOFloat* Fruiz_data = get_data_vectorf(scaling->Fruiz); + QOCOFloat* bdata = get_data_vectorf(data->b); + QOCOFloat* hdata = get_data_vectorf(data->h); + + QOCOFloat* ones_n = (QOCOFloat*)qoco_malloc(sizeof(QOCOFloat) * data->n); + for (QOCOInt j = 0; j < data->n; ++j) { + ones_n[j] = 1.0; + } + + QOCOFloat tol = 1e-6; + + // Set E[i] = 1/||row_i(A)||_inf and scale rows of A by E[i]. + if (get_nnz(data->A) > 0) { + row_inf_norm_matrix(data->A, Eruiz_data); + for (QOCOInt i = 0; i < data->p; ++i) { + Eruiz_data[i] = (Eruiz_data[i] > tol) ? (1.0 / Eruiz_data[i]) : 1.0 / tol; + } + row_col_scale_matrix(data->A, Eruiz_data, ones_n); + row_col_scale_matrix(data->At, ones_n, Eruiz_data); + } + + // Set F[i] = 1/||row_i(G)||_inf for inequality rows and scale those rows of + // G. + if (get_nnz(data->G) > 0 && data->l > 0) { + row_inf_norm_matrix(data->G, Fruiz_data); + for (QOCOInt i = 0; i < data->l; ++i) { + Fruiz_data[i] = (Fruiz_data[i] > tol) ? (1.0 / Fruiz_data[i]) : 1.0 / tol; + } + for (QOCOInt i = data->l; i < data->m; ++i) { + Fruiz_data[i] = 1.0; + } + row_col_scale_matrix(data->G, Fruiz_data, ones_n); + row_col_scale_matrix(data->Gt, ones_n, Fruiz_data); + } + + qoco_free(ones_n); + + // Scale b and h. + ew_product(bdata, Eruiz_data, bdata, data->p); + ew_product(hdata, Fruiz_data, hdata, data->m); + + // Compute inverses. + reciprocal_vectorf(scaling->Eruiz, scaling->Einvruiz); + reciprocal_vectorf(scaling->Fruiz, scaling->Finvruiz); + + set_cpu_mode(0); + + // Sync updated scaling vectors and scaled data to device (CUDA backend). + if (data->p > 0) { + sync_matrix_to_device(data->A); + sync_matrix_to_device(data->At); + } + if (data->m > 0) { + sync_matrix_to_device(data->G); + sync_matrix_to_device(data->Gt); + } + sync_vector_to_device(scaling->Eruiz); + sync_vector_to_device(scaling->Fruiz); + sync_vector_to_device(scaling->Einvruiz); + sync_vector_to_device(scaling->Finvruiz); +} + void unscale_variables(QOCOWorkspace* work) { ew_product(get_data_vectorf(work->x), get_data_vectorf(work->scaling->Druiz), diff --git a/src/qoco_api.c b/src/qoco_api.c index 10621206..ab205080 100644 --- a/src/qoco_api.c +++ b/src/qoco_api.c @@ -93,7 +93,7 @@ QOCOInt qoco_setup(QOCOSolver* solver, QOCOInt n, QOCOInt m, QOCOInt p, // Compute scaling statistics before equilibration and regularization. compute_scaling_statistics(data); - ruiz_equilibration(data, work->scaling, solver->settings->ruiz_iters); + row_normalization(data, work->scaling); // Regularize P. set_cpu_mode(1); @@ -405,7 +405,7 @@ void qoco_update_matrix_data(QOCOSolver* solver, QOCOFloat* Pxnew, compute_scaling_statistics(data); // Equilibrate new matrix data. - ruiz_equilibration(data, scaling, solver->settings->ruiz_iters); + row_normalization(data, scaling); // Regularize P. unregularize(Pcsc, -solver->settings->kkt_static_reg_P);