From fe04b7f74963527814d992b556ba250f7d607cec Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Sat, 11 Apr 2026 06:28:18 -0400 Subject: [PATCH] random_positive_definite_matrix.m: add condmax parameter We are getting some rare test failures due to these random positive-definite matrices having huge condition numbers. In the PCGM case... * cond(A)=52879791.988025, diff=0.014173772186 * cond(A)=4545512.611295, diff=0.000072462950 * cond(A)=5394929.793492, diff=0.000019629413 * etc. To avoid that problem consistently, we add a "condmax" parameter to the random matrix generator. The default (1000) is conservative to ensure that the tests pass even when we are not preconditioning. --- random_positive_definite_matrix.m | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/random_positive_definite_matrix.m b/random_positive_definite_matrix.m index 164cb33..f7dd34f 100644 --- a/random_positive_definite_matrix.m +++ b/random_positive_definite_matrix.m @@ -1,4 +1,4 @@ -function A = random_positive_definite_matrix(integerN, max_entry = realmax) +function A = random_positive_definite_matrix(integerN, max_entry = realmax, condmax=1000) % % Generate a random, symmetric positive-definite (SPD) matrix. % @@ -15,6 +15,10 @@ function A = random_positive_definite_matrix(integerN, max_entry = realmax) % % - ``max_entry`` -- (optional) Upper bound on the entries. % + % - ``condmax`` -- (optional, default=1000) Upper bound on + % the condition number; the matrix will be regenerated until + % its condition number is less or equal to ``condmax``. + % % OUTPUT: % % - ``A`` -- A symmetric, positive definite matrix. @@ -23,8 +27,11 @@ function A = random_positive_definite_matrix(integerN, max_entry = realmax) % provides unifrnd() pkg load statistics; - U = random_orthogonal_matrix(integerN); - d = unifrnd(eps, max_entry, 1, integerN); - D = diag(d); - A = U*D*U'; + A = zeros(integerN); % cond = Inf + while (cond(A) >= condmax) + U = random_orthogonal_matrix(integerN); + d = unifrnd(eps, max_entry, 1, integerN); + D = diag(d); + A = U*D*U'; + end end -- 2.51.0