]> gitweb.michael.orlitzky.com - octave.d.git/commitdiff
random_positive_definite_matrix.m: add condmax parameter
authorMichael Orlitzky <michael@orlitzky.com>
Sat, 11 Apr 2026 10:28:18 +0000 (06:28 -0400)
committerMichael Orlitzky <michael@orlitzky.com>
Sun, 12 Apr 2026 01:06:45 +0000 (21:06 -0400)
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

index 164cb33853791b7af132563274b032db41b3ba9b..f7dd34f70eb746aee2ae5b458811ee1cd7e42a6d 100644 (file)
@@ -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