From f9d6678d554de393071e8e5a5ad0972d500eaefc Mon Sep 17 00:00:00 2001 From: Ruifeng Zheng Date: Thu, 17 Sep 2026 14:08:33 +0000 Subject: [PATCH 1/4] [SPARK-59618][ML] Prevent concurrent F2J LAPACK initialization --- .../spark/ml/linalg/LAPACKInitializer.scala | 36 ++++++++++++++++++ .../distribution/MultivariateGaussian.scala | 5 ++- .../MultivariateGaussianSuite.scala | 38 +++++++++++++++++++ .../distribution/MultivariateGaussian.scala | 4 +- 4 files changed, 80 insertions(+), 3 deletions(-) create mode 100644 mllib-local/src/main/scala/org/apache/spark/ml/linalg/LAPACKInitializer.scala diff --git a/mllib-local/src/main/scala/org/apache/spark/ml/linalg/LAPACKInitializer.scala b/mllib-local/src/main/scala/org/apache/spark/ml/linalg/LAPACKInitializer.scala new file mode 100644 index 0000000000000..5fbd0bac5fc94 --- /dev/null +++ b/mllib-local/src/main/scala/org/apache/spark/ml/linalg/LAPACKInitializer.scala @@ -0,0 +1,36 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +package org.apache.spark.ml.linalg + +import breeze.linalg.{eigSym, DenseMatrix => BDM} + +private[spark] object LAPACKInitializer { + + // F2J's machine-parameter discovery uses unsynchronized mutable static state. A concurrent + // first invocation can observe partially initialized values and make some LAPACK routines loop + // forever. JVM class initialization serializes this warm-up before real calls run concurrently. + // It runs once per Spark classloader in each executor JVM: the first thread performs the warm-up + // while concurrent threads wait for class initialization to finish. Later initialize() calls are + // cheap reads of the completed Unit field. In local mode, the driver and executor share the JVM, + // so the warm-up also runs only once there. + private val initialized: Unit = { + eigSym(BDM.eye[Double](2)) + } + + def initialize(): Unit = initialized +} diff --git a/mllib-local/src/main/scala/org/apache/spark/ml/stat/distribution/MultivariateGaussian.scala b/mllib-local/src/main/scala/org/apache/spark/ml/stat/distribution/MultivariateGaussian.scala index 42746b5727029..329fb24b66cba 100644 --- a/mllib-local/src/main/scala/org/apache/spark/ml/stat/distribution/MultivariateGaussian.scala +++ b/mllib-local/src/main/scala/org/apache/spark/ml/stat/distribution/MultivariateGaussian.scala @@ -21,7 +21,7 @@ import breeze.linalg.{diag, eigSym, max, DenseMatrix => BDM, DenseVector => BDV} import org.apache.spark.annotation.{DeveloperApi, Since} import org.apache.spark.ml.impl.Utils -import org.apache.spark.ml.linalg._ +import org.apache.spark.ml.linalg.{LAPACKInitializer, _} /** @@ -108,8 +108,9 @@ class MultivariateGaussian @Since("2.0.0") ( * pseudo-determinant and the pseudo-inverse (Moore-Penrose). Singular values are considered * to be non-zero only if they exceed a tolerance based on machine precision, matrix size, and * relation to the maximum singular value (same tolerance used by, e.g., Octave). - */ + */ private def calculateCovarianceConstants: (BDM[Double], Double) = { + LAPACKInitializer.initialize() val eigSym.EigSym(d, u) = eigSym(cov.asBreeze.toDenseMatrix) // sigma = u * diag(d) * u.t // For numerical stability, values are considered to be non-zero only if they exceed tol. diff --git a/mllib-local/src/test/scala/org/apache/spark/ml/stat/distribution/MultivariateGaussianSuite.scala b/mllib-local/src/test/scala/org/apache/spark/ml/stat/distribution/MultivariateGaussianSuite.scala index f2ecff1cc58bd..d9ce968c9873d 100644 --- a/mllib-local/src/test/scala/org/apache/spark/ml/stat/distribution/MultivariateGaussianSuite.scala +++ b/mllib-local/src/test/scala/org/apache/spark/ml/stat/distribution/MultivariateGaussianSuite.scala @@ -17,6 +17,11 @@ package org.apache.spark.ml.stat.distribution +import java.util.concurrent.{CountDownLatch, Executors, TimeUnit} + +import scala.concurrent.{Await, ExecutionContext, Future} +import scala.concurrent.duration._ + import org.apache.spark.ml.SparkMLFunSuite import org.apache.spark.ml.linalg.{Matrices, Vectors} import org.apache.spark.ml.util.TestingUtils._ @@ -24,6 +29,39 @@ import org.apache.spark.ml.util.TestingUtils._ class MultivariateGaussianSuite extends SparkMLFunSuite { + test("SPARK-59618: multivariate Gaussian can be initialized concurrently") { + val numThreads = 16 + val ready = new CountDownLatch(numThreads) + val start = new CountDownLatch(1) + val executor = Executors.newFixedThreadPool(numThreads, (runnable: Runnable) => { + val thread = new Thread(runnable, "multivariate-gaussian-test") + thread.setDaemon(true) + thread + }) + implicit val executionContext: ExecutionContext = ExecutionContext.fromExecutorService(executor) + + try { + val results = (0 until numThreads).map { _ => + Future { + ready.countDown() + start.await() + val dist = new MultivariateGaussian( + Vectors.dense(0.0, 0.0), + Matrices.dense(2, 2, Array(1.0, 0.0, 0.0, 1.0))) + dist.pdf(Vectors.dense(0.0, 0.0)) + } + } + assert(ready.await(10, TimeUnit.SECONDS)) + start.countDown() + + Await.result(Future.sequence(results), 30.seconds).foreach { result => + assert(result ~== 0.15915 absTol 1E-5) + } + } finally { + executor.shutdownNow() + } + } + test("univariate") { val x1 = Vectors.dense(0.0) val x2 = Vectors.dense(1.5) diff --git a/mllib/src/main/scala/org/apache/spark/mllib/stat/distribution/MultivariateGaussian.scala b/mllib/src/main/scala/org/apache/spark/mllib/stat/distribution/MultivariateGaussian.scala index 310eb8639b3c7..eee7ee9d7bf0e 100644 --- a/mllib/src/main/scala/org/apache/spark/mllib/stat/distribution/MultivariateGaussian.scala +++ b/mllib/src/main/scala/org/apache/spark/mllib/stat/distribution/MultivariateGaussian.scala @@ -20,6 +20,7 @@ package org.apache.spark.mllib.stat.distribution import breeze.linalg.{diag, eigSym, max, DenseMatrix => DBM, DenseVector => DBV, Vector => BV} import org.apache.spark.annotation.Since +import org.apache.spark.ml.linalg.LAPACKInitializer import org.apache.spark.mllib.linalg.{Matrices, Matrix, Vector, Vectors} import org.apache.spark.mllib.util.MLUtils @@ -117,8 +118,9 @@ class MultivariateGaussian @Since("1.3.0") ( * pseudo-determinant and the pseudo-inverse (Moore-Penrose). Singular values are considered * to be non-zero only if they exceed a tolerance based on machine precision, matrix size, and * relation to the maximum singular value (same tolerance used by, e.g., Octave). - */ + */ private def calculateCovarianceConstants: (DBM[Double], Double) = { + LAPACKInitializer.initialize() val eigSym.EigSym(d, u) = eigSym(sigma.asBreeze.toDenseMatrix) // sigma = u * diag(d) * u.t // For numerical stability, values are considered to be non-zero only if they exceed tol. From 449c6927af632aa00d7a7bbd7ec520e0bc8e6ae2 Mon Sep 17 00:00:00 2001 From: Ruifeng Zheng Date: Thu, 17 Sep 2026 14:27:48 +0000 Subject: [PATCH 2/4] [SPARK-59618][ML] Refine concurrent initialization test --- .../distribution/MultivariateGaussian.scala | 2 +- .../MultivariateGaussianSuite.scala | 28 +++++++++---------- 2 files changed, 14 insertions(+), 16 deletions(-) diff --git a/mllib-local/src/main/scala/org/apache/spark/ml/stat/distribution/MultivariateGaussian.scala b/mllib-local/src/main/scala/org/apache/spark/ml/stat/distribution/MultivariateGaussian.scala index 329fb24b66cba..e74b261d0974b 100644 --- a/mllib-local/src/main/scala/org/apache/spark/ml/stat/distribution/MultivariateGaussian.scala +++ b/mllib-local/src/main/scala/org/apache/spark/ml/stat/distribution/MultivariateGaussian.scala @@ -21,7 +21,7 @@ import breeze.linalg.{diag, eigSym, max, DenseMatrix => BDM, DenseVector => BDV} import org.apache.spark.annotation.{DeveloperApi, Since} import org.apache.spark.ml.impl.Utils -import org.apache.spark.ml.linalg.{LAPACKInitializer, _} +import org.apache.spark.ml.linalg._ /** diff --git a/mllib-local/src/test/scala/org/apache/spark/ml/stat/distribution/MultivariateGaussianSuite.scala b/mllib-local/src/test/scala/org/apache/spark/ml/stat/distribution/MultivariateGaussianSuite.scala index d9ce968c9873d..835e9727176eb 100644 --- a/mllib-local/src/test/scala/org/apache/spark/ml/stat/distribution/MultivariateGaussianSuite.scala +++ b/mllib-local/src/test/scala/org/apache/spark/ml/stat/distribution/MultivariateGaussianSuite.scala @@ -17,10 +17,7 @@ package org.apache.spark.ml.stat.distribution -import java.util.concurrent.{CountDownLatch, Executors, TimeUnit} - -import scala.concurrent.{Await, ExecutionContext, Future} -import scala.concurrent.duration._ +import java.util.concurrent.{Callable, CountDownLatch, Executors, TimeUnit} import org.apache.spark.ml.SparkMLFunSuite import org.apache.spark.ml.linalg.{Matrices, Vectors} @@ -38,24 +35,25 @@ class MultivariateGaussianSuite extends SparkMLFunSuite { thread.setDaemon(true) thread }) - implicit val executionContext: ExecutionContext = ExecutionContext.fromExecutorService(executor) try { val results = (0 until numThreads).map { _ => - Future { - ready.countDown() - start.await() - val dist = new MultivariateGaussian( - Vectors.dense(0.0, 0.0), - Matrices.dense(2, 2, Array(1.0, 0.0, 0.0, 1.0))) - dist.pdf(Vectors.dense(0.0, 0.0)) - } + executor.submit(new Callable[Double] { + override def call(): Double = { + ready.countDown() + start.await() + val dist = new MultivariateGaussian( + Vectors.dense(0.0, 0.0), + Matrices.dense(2, 2, Array(1.0, 0.0, 0.0, 1.0))) + dist.pdf(Vectors.dense(0.0, 0.0)) + } + }) } assert(ready.await(10, TimeUnit.SECONDS)) start.countDown() - Await.result(Future.sequence(results), 30.seconds).foreach { result => - assert(result ~== 0.15915 absTol 1E-5) + results.foreach { result => + assert(result.get(30, TimeUnit.SECONDS) ~== 0.15915 absTol 1E-5) } } finally { executor.shutdownNow() From 007857c08ab67142beb93e9f553a4bc14a58e7f4 Mon Sep 17 00:00:00 2001 From: Ruifeng Zheng Date: Thu, 17 Sep 2026 14:36:09 +0000 Subject: [PATCH 3/4] [SPARK-59618][ML] Clarify LAPACK initialization behavior --- .../org/apache/spark/ml/linalg/LAPACKInitializer.scala | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/mllib-local/src/main/scala/org/apache/spark/ml/linalg/LAPACKInitializer.scala b/mllib-local/src/main/scala/org/apache/spark/ml/linalg/LAPACKInitializer.scala index 5fbd0bac5fc94..6d4e61ec08469 100644 --- a/mllib-local/src/main/scala/org/apache/spark/ml/linalg/LAPACKInitializer.scala +++ b/mllib-local/src/main/scala/org/apache/spark/ml/linalg/LAPACKInitializer.scala @@ -28,6 +28,12 @@ private[spark] object LAPACKInitializer { // while concurrent threads wait for class initialization to finish. Later initialize() calls are // cheap reads of the completed Unit field. In local mode, the driver and executor share the JVM, // so the warm-up also runs only once there. + // + // The warm-up uses Breeze's selected LAPACK backend. It is required only for F2J; a native + // provider such as OpenBLAS or MKL only incurs this one-time 2-by-2 decomposition. Later calls + // remain concurrent. A 2-by-2 matrix is used because DSYEV returns early for a 1-by-1 matrix + // without initializing floating-point limits such as epsilon, the safe minimum, and the maximum + // finite value. private val initialized: Unit = { eigSym(BDM.eye[Double](2)) } From ffd09d1313b1a71b84fd2852aca7691a1da0cc9a Mon Sep 17 00:00:00 2001 From: Ruifeng Zheng Date: Thu, 17 Sep 2026 14:41:12 +0000 Subject: [PATCH 4/4] [SPARK-59618][ML] Explain F2J machine parameters --- .../spark/ml/linalg/LAPACKInitializer.scala | 20 +++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/mllib-local/src/main/scala/org/apache/spark/ml/linalg/LAPACKInitializer.scala b/mllib-local/src/main/scala/org/apache/spark/ml/linalg/LAPACKInitializer.scala index 6d4e61ec08469..913048b5c567b 100644 --- a/mllib-local/src/main/scala/org/apache/spark/ml/linalg/LAPACKInitializer.scala +++ b/mllib-local/src/main/scala/org/apache/spark/ml/linalg/LAPACKInitializer.scala @@ -21,19 +21,23 @@ import breeze.linalg.{eigSym, DenseMatrix => BDM} private[spark] object LAPACKInitializer { - // F2J's machine-parameter discovery uses unsynchronized mutable static state. A concurrent - // first invocation can observe partially initialized values and make some LAPACK routines loop - // forever. JVM class initialization serializes this warm-up before real calls run concurrently. + // Breeze's eigSym delegates to the selected LAPACK backend. When LAPACK selects its F2J fallback, + // F2J's Dlamc2 discovers and caches machine parameters: the floating-point radix and precision, + // rounding behavior, epsilon, minimum safe positive value and exponent, and maximum finite value + // and exponent. Dlamc2 stores them in unsynchronized mutable static fields and sets first to + // false before all fields are initialized. A concurrent first invocation can therefore observe + // partial values and make some LAPACK routines loop forever. + // + // JVM class initialization serializes this warm-up before real calls run concurrently. // It runs once per Spark classloader in each executor JVM: the first thread performs the warm-up // while concurrent threads wait for class initialization to finish. Later initialize() calls are // cheap reads of the completed Unit field. In local mode, the driver and executor share the JVM, // so the warm-up also runs only once there. // - // The warm-up uses Breeze's selected LAPACK backend. It is required only for F2J; a native - // provider such as OpenBLAS or MKL only incurs this one-time 2-by-2 decomposition. Later calls - // remain concurrent. A 2-by-2 matrix is used because DSYEV returns early for a 1-by-1 matrix - // without initializing floating-point limits such as epsilon, the safe minimum, and the maximum - // finite value. + // The warm-up follows the same Breeze -> LAPACK -> F2J path as the real operation when F2J is + // selected. A native provider such as OpenBLAS or MKL only incurs this one-time 2-by-2 + // decomposition, and later calls remain concurrent. A 2-by-2 matrix is used because DSYEV returns + // early for a 1-by-1 matrix without requesting the machine parameters. private val initialized: Unit = { eigSym(BDM.eye[Double](2)) }