WIP -- added placeholder file for the nearPD function #800
Conversation
| * @author ewalsh | ||
| **/ | ||
|
|
||
| class NearPD(matA: DenseMatrix[Double], corrBool: Boolean = false, keepDiag: Boolean = false, |
There was a problem hiding this comment.
Is there a reason this needs to be a class? AFAICT you're only doing it to make it top level? Could you maybe put it in the linalg package object as a def? or make it an object with an apply method?
There was a problem hiding this comment.
(I don't think there's a reason to make this a "UFunc" fwiw)
There was a problem hiding this comment.
First, thank you for your comments. I found all of them very helpful, and they are very much appreciated. It may take me a hot minute but I will address everything ASAP.
There was a problem hiding this comment.
I have addressed or commented on each of these points. Please feel free to add any other comments, I will defer to you to revolve each based on your approval.
I have changed to an object with an apply method
| * applied to the result of the Higham algorithm. | ||
| * @param doSym Boolean indicating if X <- (X + t(X))/2 should be done, | ||
| * after X <- tcrossprod(Qd, Q); some doubt if this is necessary. | ||
| * @param doDykstra Boolean indicating if Dykstra's correction should be used; |
There was a problem hiding this comment.
i'd prefer there be an algo: NearPD.Algorithm or something, where NearPD.Algorithm is a sealed trait or an enum
There was a problem hiding this comment.
I have added sealed traits for both the algorithm and choice of convergence norm type. This makes more sense and I hope is easier to read.
| //--- start algorithm --- | ||
| def generate: DenseMatrix[Double] = { | ||
| // choose Dykstra or direct | ||
| val matX = chooseAlgo(symA, doDykstra) |
There was a problem hiding this comment.
I don't like this. how about inlining?
val matX = if(algo == Dystrka) doDystrkaFunc( ) else doDirect()
There was a problem hiding this comment.
I have changed style to inline in each opportunity that I found
| } | ||
| } | ||
| // keepDiag function | ||
| def keepDiagFunc(mat: DenseMatrix[Double], bool: Boolean): DenseMatrix[Double] = { |
There was a problem hiding this comment.
as with chooseAlgo i'd rather you just inline the boolean gating of this method and this method only does the keepDiag part
| } | ||
| } | ||
| // doSym function | ||
| def doSymFunc(mat: DenseMatrix[Double], bool: Boolean): DenseMatrix[Double] = { |
| val matQ = fliplr(eigVecVals.eigenvectors) | ||
| val vecd = reverse(eigVecVals.eigenvalues) | ||
| val p = vecd.map(x => x > eigTol) | ||
| val pcheck = breeze.linalg.sum(p.map(x => bool2Dbl(x))) |
| val pcheck = breeze.linalg.sum(p.map(x => bool2Dbl(x))) | ||
| if(pcheck < 1.0){ | ||
| println("Matrix seems negative semi-definite") | ||
| break |
There was a problem hiding this comment.
you have to wrap the loop you're breaking out of in breakable in Scala. It's dumb
| matX | ||
| } | ||
|
|
||
| def doDirectFixedPt(mat: DenseMatrix[Double]): DenseMatrix[Double] = { |
There was a problem hiding this comment.
this function shares a ton of code with the one above. Is it possible to extra out more. I know I had you inline some stuff, but I think it would be cleaner to extract out the pcheck stuff, and the eigVecVals etc stuff, or possibly to make it so that the main loop takes a function as input that does the differences between the two algorithms, though I think that might be worse.
There was a problem hiding this comment.
This one is the only one I haven't truly addressed. The rationale is this, it would be possible to extract out a lot of the code here. However, when running the function it would require running an if statement within the loop to check which algorithm to run.
Please advise. I'd be happy to change, but I just wanted to point out the original rationale. I think in large cases where lots of iterations are needed this repetitive checking could be slower than necessary.
| // ## A longer example, extended from Jens' original, | ||
| // ## showing the effects of some of the options: | ||
|
|
||
| // val matA = DenseMatrix(Array(Array(1, 0.477, 0.644, 0.478, 0.651, 0.826), |
There was a problem hiding this comment.
it would be good to make this into a test, and maybe add some randomized tests about properties this function should have, e.g.: for any matrix you give it (that one might want to give it), the output is a PD matrix or it errors, and it doesn't error too much, or something like that.
There was a problem hiding this comment.
Done and I also added another simple test
| } | ||
|
|
||
| // bool to double | ||
| def bool2Dbl(bool: Boolean): Double = { |
There was a problem hiding this comment.
i don't think you actually need this here, but if you do, there's a function in breeze called I for indicator function that does exactly this, except it works on vectors
| converged = (conv < convTol) | ||
| } | ||
| if(!converged){ | ||
| println(s"Did not converge in ${iter} iterations") |
| //--- start algorithm --- | ||
| def generate: DenseMatrix[Double] = { | ||
| // choose Dykstra or direct | ||
| @inline val matX = if(algo == Dykstra) { doDykstraFunc(symA) } else { doDirectFixedPt(symA) } |
| val n = mat.cols - 1 | ||
| val eps = posdTol * breeze.numerics.abs(vecd(0)) | ||
| if(vecd(n) < eps){ | ||
| vecd = vecd.map(x => { |
| } | ||
|
|
||
| // norm type function | ||
| def normTypeFunc(matX: DenseMatrix[Double], matY: DenseMatrix[Double], normType: ConvergenceNormType): Double = { |
There was a problem hiding this comment.
i should really fix norm for matrices...
| // break if negative semi-definite | ||
| def negSemiDefCheckBreak(pcheck: Boolean): Unit = { | ||
| if(!pcheck){ | ||
| println("Matrix seems negative semi-definite") |
Adding a function to find the nearest positive definite matrix via the Higham (2002) algorithm.