Skip to content

Commit 0140a4c

Browse files
jajhalljajhall
committed
Got Power method right, but no options not correct
1 parent 1b836fa commit 0140a4c

File tree

4 files changed

+137
-20
lines changed

4 files changed

+137
-20
lines changed

‎check/TestPdlp.cpp‎

Lines changed: 12 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -331,17 +331,26 @@ TEST_CASE("pdlp-restart-add-row", "[pdlp]") {
331331
}
332332

333333
TEST_CASE("hi-pdlp", "[pdlp]") {
334-
std::string model = "adlittle";//"avgas";
334+
std::string model = "avgas";//"adlittle";//"afiro";//
335335
std::string model_file =
336336
std::string(HIGHS_DIR) + "/check/instances/" + model + ".mps";
337337
Highs h;
338338
// h.setOptionValue("output_flag", dev_run);
339339
REQUIRE(h.readModel(model_file) == HighsStatus::kOk);
340340
h.setOptionValue("solver", kHiPdlpString);
341341
h.setOptionValue("kkt_tolerance", kkt_tolerance);
342-
h.setOptionValue("pdlp_iteration_limit", 10000);
343-
h.setOptionValue("log_dev_level", kHighsLogDevLevelVerbose);
342+
h.setOptionValue("pdlp_scaling_mode", 0);
343+
h.setOptionValue("pdlp_restart_strategy", 0);
344+
h.setOptionValue("pdlp_step_size_strategy", 0);
345+
// h.setOptionValue("pdlp_iteration_limit", 10000);
346+
// h.setOptionValue("log_dev_level", kHighsLogDevLevelVerbose);
344347
HighsStatus run_status = h.run();
345348
// REQUIRE(run_status == HighsStatus::kOk);
346349
// REQUIRE(h.getModelStatus() == HighsModelStatus::kOptimal);
350+
const bool cupdlp_test = false;
351+
if (cupdlp_test) {
352+
h.setOptionValue("solver", kCuPdlpString);
353+
h.setOptionValue("pdlp_features_off", kPdlpAllFeaturesOff);
354+
run_status = h.run();
355+
}
347356
}

‎highs/pdlp/cupdlp/cupdlp_step.c‎

Lines changed: 34 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -77,6 +77,8 @@ cupdlp_retcode PDHG_Power_Method(CUPDLPwork *work, cupdlp_float *lambda) {
7777
if (work->settings->nLogLevel>0)
7878
cupdlp_printf("Power Method:\n");
7979

80+
// work->buffer is a CUPDLPvec*, but a pointer to its values is
81+
// needed for some linalg calls
8082
cupdlp_float *q = work->buffer->data;
8183

8284
cupdlp_initvec(q, 1.0, lp->nRows);
@@ -86,9 +88,24 @@ cupdlp_retcode PDHG_Power_Method(CUPDLPwork *work, cupdlp_float *lambda) {
8688
CUPDLPvec *aty = iterates->aty[iter % 2];
8789

8890
double res = 0.0;
91+
double previous_lambda = 0.0;
92+
// q is stored in work->buffer->data
93+
//
94+
// z is stored in ax->data
95+
//
96+
// y is stored in aty->data
97+
//
98+
int log_iters = work->settings->nLogLevel > 0;
99+
log_iters = 1;
100+
101+
if (log_iters)
102+
cupdlp_printf("It lambda dl_lambda residual\n");
89103
for (cupdlp_int iter = 0; iter < 20; ++iter) {
90104
// z = A*A'*q
105+
//
106+
// as y = A'q...
91107
ATy(work, aty, work->buffer);
108+
// ... then z = Ay
92109
Ax(work, ax, aty);
93110

94111
// q = z / norm(z)
@@ -97,17 +114,30 @@ cupdlp_retcode PDHG_Power_Method(CUPDLPwork *work, cupdlp_float *lambda) {
97114
cupdlp_twoNorm(work, lp->nRows, q, &qNorm);
98115
cupdlp_scaleVector(work, 1.0 / qNorm, q, lp->nRows);
99116

117+
// Now compute the Rayleigh quotient of q which, since q'q=1 is
118+
//
119+
// lambda = q'AA'q = w^Tw = ||w||_2
120+
//
121+
// where w = A'q
122+
//
123+
// aty is no longer needed, so w is stored in aty->data
100124
ATy(work, aty, work->buffer);
101-
102125
cupdlp_twoNormSquared(work, lp->nCols, aty->data, lambda);
103126

104127
cupdlp_float alpha = -(*lambda);
128+
// Now compute the residual between z = A*A'*q (old q) and
129+
// lambda.q (new q) to just to log progress
130+
//
131+
// z := z - lambda.q
105132
cupdlp_axpy(work, lp->nRows, &alpha, q, ax->data);
106133

107134
cupdlp_twoNormSquared(work, lp->nCols, ax->data, &res);
108135

109-
if (work->settings->nLogLevel>0)
110-
cupdlp_printf("% d %e %.3f\n", iter, *lambda, res);
136+
double dl_lambda = fabs(*lambda - previous_lambda);
137+
previous_lambda = *lambda;
138+
139+
if (log_iters)
140+
cupdlp_printf("%2d %12.6g %11.4g %11.4g\n", iter, *lambda, dl_lambda, res);
111141
}
112142

113143
exit_cleanup:
@@ -424,4 +454,4 @@ cupdlp_retcode PDHG_Update_Iterate(CUPDLPwork *pdhg) {
424454

425455
exit_cleanup:
426456
return RETCODE_OK;
427-
}
457+
}

‎highs/pdlp/hipdlp/logger.hpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -39,7 +39,7 @@ class Logger {
3939
public:
4040
Logger(LogLevel level = LogLevel::kInfo);
4141
void set_log_file(const std::string& filename);
42-
42+
LogLevel getLogLevel() const { return console_level_; }
4343
// Logging methods for different levels
4444
void info(const std::string& message);
4545
void verbose(const std::string& message);

‎highs/pdlp/hipdlp/pdhg.cc‎

Lines changed: 90 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -792,40 +792,118 @@ HighsStatus PDLPSolver::PowerMethod(HighsLp &lp, double& op_norm_sq) {
792792
const double tol = 1e-6;
793793

794794
// Initialize a random vector x
795-
std::vector<double> x(lp.num_col_);
795+
std::vector<double> x_vec(lp.num_col_);
796796
std::random_device rd;
797797
std::mt19937 engine_fixed_seed(12345); //gen(rd());
798798
std::uniform_real_distribution<> dis(-1.0, 1.0);
799799
for (HighsInt i = 0; i < lp.num_col_; ++i) {
800-
x[i] = dis(engine_fixed_seed);
800+
x_vec[i] = dis(engine_fixed_seed);
801801
}
802-
linalg::normalize(x); // Assumes a normalize function in linalg
802+
linalg::normalize(x_vec); // Assumes a normalize function in linalg
803803

804+
const HighsInt kYanyuPowerMethod = 0;
805+
const HighsInt kATAPowerMethod = 1;
806+
const HighsInt kCuPdlpAATPowerMethod = 2;
807+
808+
const HighsInt power_method = kCuPdlpAATPowerMethod;
804809
// Allocate memory for matrix-vector products
805-
std::vector<double> Ax_vec(lp.num_row_);
806-
std::vector<double> ATAx_vec(lp.num_col_);
810+
std::vector<double> y_vec;
811+
std::vector<double> z_vec;
812+
if (power_method < kCuPdlpAATPowerMethod) {
813+
y_vec.resize(lp.num_row_);
814+
z_vec.resize(lp.num_col_);
815+
} else {
816+
y_vec.resize(lp.num_col_);
817+
z_vec.resize(lp.num_row_);
818+
}
807819

808820
double op_norm_sq_old = 0.0;
821+
LogLevel log_level = logger_.getLogLevel();
822+
int log_iters = log_level == LogLevel::kVerbose || log_level == LogLevel::kDebug;
823+
log_iters = 1;
824+
825+
if (log_iters) printf("It lambda dl_lambda\n");
809826

827+
if (power_method == kATAPowerMethod) {
828+
x_vec.assign(lp.num_col_, 1);
829+
} else if (power_method == kCuPdlpAATPowerMethod) {
830+
x_vec.assign(lp.num_row_, 1);
831+
}
832+
const HighsSparseMatrix& matrix = lp.a_matrix_;
833+
double lambda = 0.0;
834+
double previous_lambda = lambda;
810835
for (int iter = 0; iter < max_iter; ++iter) {
836+
if (power_method == kYanyuPowerMethod) {
837+
// Original Yanyu power method
838+
//
811839
// Compute ATAx = A^T * (A * x)
812-
linalg::Ax(lp, x, Ax_vec);
813-
linalg::ATy(lp, Ax_vec, ATAx_vec); // Note: ATy computes A^T * vector
840+
linalg::Ax(lp, x_vec, y_vec);
841+
linalg::ATy(lp, y_vec, z_vec); // Note: ATy computes A^T * vector
814842

815843
// Estimate the squared operator norm (largest eigenvalue of A^T*A)
816-
op_norm_sq = linalg::dot(x, ATAx_vec); // Assumes a dot product function in linalg
844+
op_norm_sq = linalg::dot(x_vec, z_vec); // Assumes a dot product function in linalg
817845

818846
// Check for convergence
819847
if (std::abs(op_norm_sq - op_norm_sq_old) < tol * op_norm_sq) {
820848
return HighsStatus::kOk;
821849
}
850+
double dl_op_norm_sq = std::fabs(op_norm_sq - op_norm_sq_old);
822851
op_norm_sq_old = op_norm_sq;
823852

824853
// Prepare for the next iteration
825-
linalg::normalize(ATAx_vec); // Normalize the result
826-
x = ATAx_vec;
827-
}
828-
854+
linalg::normalize(z_vec); // Normalize the result
855+
x_vec = z_vec;
856+
if (log_iters) printf("%2d %12.6g %11.4g\n", iter, op_norm_sq, dl_op_norm_sq);
857+
} else {
858+
if (power_method == kATAPowerMethod) {
859+
// Yanyu power method corrected - with Rayleigh quotient
860+
//
861+
// Compute z = ATAx = A^T * (A * x)
862+
matrix.product(y_vec, x_vec);
863+
matrix.productTranspose(z_vec, y_vec);
864+
// q = z / norm(z)
865+
double z_norm = std::sqrt(linalg::dot(z_vec, z_vec));
866+
double denom = 1.0 / z_norm;
867+
for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++)
868+
z_vec[iCol] *= denom;
869+
// Now compute the Rayleigh quotient of q which, since q'q=1 is
870+
//
871+
// lambda = q'A'Aq = w^Tw = ||w||_2
872+
//
873+
// where w = Aq
874+
//
875+
// y_vec is no longer needed, so w is stored in it
876+
matrix.product(y_vec, z_vec);
877+
lambda = linalg::dot(y_vec, y_vec);
878+
} else {
879+
// cuPDLP-C power method
880+
//
881+
// Compute z = AA^Tx = A * (A^T * x)
882+
matrix.productTranspose(y_vec, x_vec);
883+
matrix.product(z_vec, y_vec);
884+
// q = z / norm(z)
885+
double z_norm = std::sqrt(linalg::dot(z_vec, z_vec));
886+
double denom = 1.0 / z_norm;
887+
for (HighsInt iRow = 0; iRow < lp.num_row_; iRow++)
888+
z_vec[iRow] *= denom;
889+
// Now compute the Rayleigh quotient of q which, since q'q=1 is
890+
//
891+
// lambda = q'AA'q = w^Tw = ||w||_2
892+
//
893+
// where w = A'q
894+
//
895+
// y_vec is no longer needed, so w is stored in it
896+
matrix.productTranspose(y_vec, z_vec);
897+
lambda = linalg::dot(y_vec, y_vec);
898+
}
899+
double dl_lambda = std::fabs(lambda - previous_lambda);
900+
previous_lambda = lambda;
901+
902+
x_vec = z_vec;
903+
if (log_iters) printf("%2d %12.6g %11.4g\n", iter, lambda, dl_lambda);
904+
}
905+
}
906+
if (power_method != kYanyuPowerMethod) op_norm_sq = lambda;
829907
// If the method did not converge within max_iter
830908
return HighsStatus::kWarning;
831909
}

0 commit comments

Comments
 (0)