New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
matrix multiply fails on transposed matrix, but succeeds on copy of transposed matrix #850
Comments
Here's a somewhat cleaner version of the test script, showing exactly what the problem is: import breeze.linalg._
object ExpressionCrash {
def main(args: Array[String]): Unit = {
val (matA, matB) = createTriggerData
matA * matB.copy // <<<<<< this succeeds
matA * matB // <<<<<< this crashes
}
def createTriggerData: (DenseMatrix[Double], DenseMatrix[Double]) = {
val yAndX = testmatrix
val y = yAndX(::,0) // column zero is y
val X = yAndX(::,1 until yAndX.cols) // remainder of columns are predictors
val Iᴛ: DenseMatrix[Double] = DenseMatrix.eye[Double](X.rows)
val Iɴ: DenseMatrix[Double] = DenseMatrix.eye[Double](X.cols)
val ιᴛ: DenseVector[Double] = DenseVector.ones[Double](X.rows)
val ιɴ: DenseVector[Double] = DenseVector.ones[Double](X.cols)
val T = X.rows.toDouble
val Jᴛ: DenseMatrix[Double] = Iᴛ - 1.0/T * ιᴛ * ιᴛ.t
val Jɴ: DenseMatrix[Double] = Iɴ - 1.0/T * ιɴ * ιɴ.t
val ȳ = ιᴛ.t * y/T
var Wxz: DenseVector[Double] = DenseVector.zeros(X.cols)
val Z = DenseVector.ones[Double](X.rows)
(Jɴ, X.t)
}
def J(T: Double): DenseMatrix[Double] = {
val Iᴛ: DenseMatrix[Double] = DenseMatrix.eye[Double](T.toInt)
val ιᴛ = DenseVector.ones[Double](T.toInt)
val Jᴛ = Iᴛ - 1.0/T * ιᴛ * ιᴛ.t
Jᴛ
}
lazy val (testrows, testcols) = (4, 3)
lazy val testmatrix = new DenseMatrix(testcols, testrows, testdata.toArray).t
lazy val testdata = Seq(
67, 33, 64,
62, 69, 78,
76, 58, 93,
57, 60, 60,
).map( _.toDouble ).toArray
} And here's the somewhat cleaner crash dump: $ jsrc/expressionCrash.sc Dec 02, 2022 4:21:06 PM dev.ludovic.netlib.blas.InstanceBuilder initializeNative WARNING: Failed to load implementation from:dev.ludovic.netlib.blas.JNIBLAS Dec 02, 2022 4:21:06 PM dev.ludovic.netlib.blas.InstanceBuilder initializeJava WARNING: Failed to load implementation from:dev.ludovic.netlib.blas.VectorBLAS java.lang.IndexOutOfBoundsException: Index 12 out of bounds for length 12 at dev.ludovic.netlib.blas.AbstractBLAS.checkIndex(AbstractBLAS.java:51) at dev.ludovic.netlib.blas.AbstractBLAS.dgemm(AbstractBLAS.java:295) at breeze.linalg.operators.DenseMatrixMultiplyOps$impl_OpMulMatrix_DMD_DMD_eq_DMD$.apply(DenseMatrixOps.expanded.scala:192) at breeze.linalg.operators.DenseMatrixMultiplyOps$impl_OpMulMatrix_DMD_DMD_eq_DMD$.apply(DenseMatrixOps.expanded.scala:186) at breeze.linalg.ImmutableNumericOps.$times(NumericOps.scala:119) at breeze.linalg.ImmutableNumericOps.$times$(NumericOps.scala:27) at breeze.linalg.DenseMatrix.$times(DenseMatrix.scala:52) at ExpressionCrash$.main(expressionCrash.sc:10) at ExpressionCrash.main(expressionCrash.sc) at java.base/jdk.internal.reflect.NativeMethodAccessorImpl.invoke0(Native Method) at java.base/jdk.internal.reflect.NativeMethodAccessorImpl.invoke(NativeMethodAccessorImpl.java:78) at java.base/jdk.internal.reflect.DelegatingMethodAccessorImpl.invoke(DelegatingMethodAccessorImpl.java:43) at java.base/java.lang.reflect.Method.invoke(Method.java:567) at dotty.tools.scripting.ScriptingDriver.compileAndRun(ScriptingDriver.scala:33) at dotty.tools.scripting.Main$.process(Main.scala:45) at dotty.tools.MainGenericRunner$.run$1(MainGenericRunner.scala:250) at dotty.tools.MainGenericRunner$.process(MainGenericRunner.scala:270) at dotty.tools.MainGenericRunner$.main(MainGenericRunner.scala:281) at dotty.tools.MainGenericRunner.main(MainGenericRunner.scala) BTW, it seems possible that this is a bug in dev.ludovic.netlib.blas.AbstractBLAS.checkIndex rather than breeze. On line 295, the method is called with index==12, which is clearly beyond the end of the array, but the actual expression would seem to only require 11. I opened an issue here: netlib issue 19 |
This seems like a bug to me, hopefully it's not a misunderstanding on my part.
The expression on line 35 executes as expected, but the seemingly identical expression on line 38 causes a runtime crash. The main difference between the expression
X.t
and(X.t).copy
AFAICT is the tranposed matrix has an offset of 1 column into a larger underlying array, due to the way it was created.The behavior is the same with both scala 2.13.10 and scala 3.2.1.
I have tested in the following environments:
The code:
Output with stack dump:
The text was updated successfully, but these errors were encountered: