Skip to content
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

WIP -- added placeholder file for the nearPD function #800

Open
wants to merge 5 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from 4 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -7,3 +7,4 @@ project
.*.swp
.idea*
.idea
.classpath
276 changes: 276 additions & 0 deletions math/src/main/scala/breeze/linalg/functions/nearpd.scala
Original file line number Diff line number Diff line change
@@ -0,0 +1,276 @@
// WIP file holder for a find nearest positive definite matrix function

package breeze.linalg.functions

import breeze.linalg.DenseMatrix
import breeze.linalg.DenseVector
import breeze.linalg.eigSym
import breeze.linalg.fliplr
import breeze.linalg.reverse
import breeze.linalg.diag
import breeze.linalg.*
import breeze.linalg.Axis
import util.control.Breaks.break

/** Generate the nearest (in Frobenius norm) positive definite matrix
* using the algorithm decribed in (Higham 2002)
* Inspired by the NearPD function from the Matrix package
* from Martin Maechler
* @param matA DenseMatrix of doubles n * n approximately positive definite matrix,
* typically an approximation to a correlation or covariance matrix.
* If x is not symmetric (and ensureSymmetry is not false), symmpart(x) is used.
* @param corrBool Boolean indicating if the matrix should be a correlation matrix.
* @param keepDiag Boolean generalizing corr: if true, the resulting matrix
* should have the same diagonal (diag(x)) as the input matrix.
* @param do2eigen Boolean indicating if a posdefify() eigen step should be
* 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;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i'd prefer there be an algo: NearPD.Algorithm or something, where NearPD.Algorithm is a sealed trait or an enum

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

* true by default. If false, the algorithm is basically the direct fixpoint
* iteration Y(k) = P_U(P_S(Y(k-1))).
* @param onlyVals Boolean if TRUE, the result is just the vector of
* eigenvalues of the approximating matrix.
* the target vector to mean 0 and unit variance, default is false
* @param eigTol defines relative positiveness of eigenvalues compared to
* largest one, λ_1. Eigenvalues λ_k are treated as if zero when λ_k / λ_1 ≤ eig.tol.
* @param convTol convergence tolerance for Higham algorithm.
* @param posdTol tolerance for enforcing positive definiteness
* (in the final posdefify step when do2eigen is TRUE).
* @param maxit maximum number of iterations allowed.
* @param convNormType convergence norm type (norm(*, type)) used for Higham
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'd rather a sealed trait here

* algorithm. The default is "I" (infinity), for reasons of speed (and back compatibility);
* using "F" is more in line with Higham's proposal.
* @param trace Boolean specifying if convergence monitoring should be traced.
*/
/**
* TODO
* migrate to breeze-y style, integrate with library, add tests
* @author ewalsh
**/

class NearPD(matA: DenseMatrix[Double], corrBool: Boolean = false, keepDiag: Boolean = false,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

(I don't think there's a reason to make this a "UFunc" fwiw)

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

do2eigen: Boolean = true, doSym: Boolean = false, doDykstra: Boolean = true,
onlyVals: Boolean = false, eigTol: Double = 1e-6, convTol: Double = 1e-7, posdTol: Double = 1e-8,
maxit: Int = 100, convNormType: String = "I", trace: Boolean = false) {
// ensure symmetry
val symA = (matA + matA.t) / 2.0
val diagX0 = diag(symA)
//--- start algorithm ---
def generate: DenseMatrix[Double] = {
// choose Dykstra or direct
val matX = chooseAlgo(symA, doDykstra)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't like this. how about inlining?

val matX = if(algo == Dystrka) doDystrkaFunc( ) else doDirect()

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have changed style to inline in each opportunity that I found

// run posdefify ?
if(do2eigen || onlyVals){
posDefify(matX)
} else {
matX
}
}

// posdefify eigen step
def posDefify(mat: DenseMatrix[Double]): DenseMatrix[Double] = {
var matX = mat
val eigVecVals = eigSym(matX)
var vecd = reverse(eigVecVals.eigenvalues)
val n = mat.cols - 1
val eps = posdTol * breeze.numerics.abs(vecd(0))
if(vecd(n) < eps){
vecd = vecd.map(x => {
if(x < eps){
eps
} else {
x
}
})
if(!onlyVals){
val matQ = fliplr(eigVecVals.eigenvectors)
val matDiag = diag(matX)
matX = matQ * matQ.t(::, *).map(_ * vecd)
val dVec = breeze.numerics.sqrt(pmax(matDiag, DenseVector(eps))/diag(matX))
matX = (dVec * dVec.t) *:* matX
}
}
if(onlyVals){
matX = DenseMatrix(vecd)
}
if(corrBool){
matX = corrDiagFunc(matX, corrBool)
}
if(keepDiag){
matX = keepDiagFunc(matX, keepDiag)
}
matX
}
// parallel max
def pmax(vec0: DenseVector[Double], vec1: DenseVector[Double]): DenseVector[Double] = {
val max0 = breeze.linalg.max(vec0)
val max1 = breeze.linalg.max(vec1)
if(max0 > max1){
vec0
} else {
vec1
}
}
// choose algo
def chooseAlgo(mat: DenseMatrix[Double], bool: Boolean): DenseMatrix[Double] = {
bool match {
case true => doDykstraFunc(mat)
case false => doDirectFixedPt(mat)
}
}

// bool to double
def bool2Dbl(bool: Boolean): Double = {
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Agreed and removed

if(bool){
1.0
} else {
0.0
}
}
// keepDiag function
def keepDiagFunc(mat: DenseMatrix[Double], bool: Boolean): DenseMatrix[Double] = {
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

as with chooseAlgo i'd rather you just inline the boolean gating of this method and this method only does the keepDiag part

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done

bool match {
case true => {
DenseMatrix.tabulate(mat.rows, mat.cols)({
case(i, j) => {
if(i == j){
diagX0(i)
} else {
mat(i, j)
}
}
})
}
case false => mat
}
}
// doSym function
def doSymFunc(mat: DenseMatrix[Double], bool: Boolean): DenseMatrix[Double] = {
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done

bool match {
case true => (mat * mat.t) / 2.0
case false => mat
}
}
// corr diag function
def corrDiagFunc(mat: DenseMatrix[Double], bool: Boolean): DenseMatrix[Double] = {
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same

bool match {
case true => {
DenseMatrix.tabulate(mat.rows, mat.cols)({
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i might be misremembering, but I think this would be simpler/quicker as

val corrDiaged = mat.copy
diag(corrDiaged) := 1.0
corrDiaged

case(i, j) => {
if(i == j){
1.0
} else {
mat(i, j)
}
}
})
}
case false => mat
}
}
// norm type function
def normTypeFunc(matX: DenseMatrix[Double], matY: DenseMatrix[Double], normType: String): Double = {
normType match {
case "F" => breeze.linalg.norm((matY - matX).toDenseVector) / breeze.linalg.norm((matY).toDenseVector)
case "1" => breeze.linalg.max(breeze.linalg.sum(breeze.numerics.abs(matY - matX), Axis._0)) / breeze.linalg.max(breeze.linalg.sum(breeze.numerics.abs(matY), Axis._0))
case "2" => breeze.linalg.max(breeze.linalg.svd(matY - matX).singularValues) / breeze.linalg.max(breeze.linalg.svd(matY).singularValues)
case _ => breeze.linalg.max(breeze.linalg.sum(breeze.numerics.abs(matY - matX), Axis._1)) / breeze.linalg.max(breeze.linalg.sum(breeze.numerics.abs(matY), Axis._1)) // Infinity Norm
}
}
// doDykstra function
def doDykstraFunc(mat: DenseMatrix[Double]): DenseMatrix[Double] = {
var iter: Int = 0
var converged: Boolean = false
var conv = Double.PositiveInfinity
val n: Int = mat.cols
var dSubS = DenseMatrix.zeros[Double](mat.rows, n)
var matX = mat
var matY = matX
var matR = matY
while(iter < maxit && !converged){
matY = matX
matR = matY - dSubS
val eigVecVals = eigSym(matR)
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)))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

val pcheck = any(vecd > eigTol)

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

incorporated directly

if(pcheck < 1.0){
println("Matrix seems negative semi-definite")
break
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

you have to wrap the loop you're breaking out of in breakable in Scala. It's dumb

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done

}
val matQtrim = matQ(::, p).toDenseMatrix
val vecdTrim = vecd(p).toDenseVector
val cpPt1 = DenseMatrix.tabulate(matQtrim.rows, matQtrim.cols)({
case(i, j) => matQtrim(i, j) * vecdTrim(j)
})
matX = cpPt1 * matQtrim.t
dSubS = matX - matR
matX = doSymFunc(matX, doSym)
matX = corrDiagFunc(matX, corrBool)
matX = keepDiagFunc(matX, keepDiag)
conv = normTypeFunc(matX, matY, convNormType)
iter += 1
if(trace){
println(s"iter ${iter} : p = ${breeze.linalg.sum(vecd)} ||Y-X|| / ||Y|| = ${conv}")
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i prefer to use logging rather than printlns in breeze

}
converged = (conv < convTol)
}
if(!converged){
println(s"Did not converge in ${iter} iterations")
}
matX
}

def doDirectFixedPt(mat: DenseMatrix[Double]): DenseMatrix[Double] = {
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

var iter = 0
var converged = false
var conv = Double.PositiveInfinity
var matX = mat
var matY = matX
while(iter < maxit && !converged){
matY = matX
val eigVecVals = eigSym(matY)
val matQ = fliplr(eigVecVals.eigenvectors)
val vecd = reverse(eigVecVals.eigenvalues)
val p = vecd.map(x => x > eigTol).toDenseVector
val pcheck = breeze.linalg.sum(p.map(x => bool2Dbl(x)))
if(pcheck < 1.0){
println("Matrix seems negative semi-definite")
break
}
val matQtrim = matQ(::, p).toDenseMatrix
val vecdTrim = vecd(p).toDenseVector
val cpPt1 = DenseMatrix.tabulate(matQtrim.rows, matQtrim.cols)({
case(i, j) => matQtrim(i, j) * vecdTrim(j)
})
matX = cpPt1 * matQtrim.t
matX = doSymFunc(matX, doSym)
matX = corrDiagFunc(matX, corrBool)
matX = keepDiagFunc(matX, keepDiag)
conv = normTypeFunc(matX, matY, convNormType)
iter += 1
if(trace){
println(s"iter ${iter} : p = ${breeze.linalg.sum(vecd)} ||Y-X|| / ||Y|| = ${conv}")
}
converged = (conv < convTol)
}
if(!converged){
println(s"Did not converge in ${iter} iterations")
}
matX
}

}

// ## 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),
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done and I also added another simple test

// Array(0.477, 1, 0.516, 0.233, 0.682, 0.75),
// Array(0.644, 0.516, 1, 0.599, 0.581, 0.742),
// Array(0.478, 0.233, 0.599, 1, 0.741, 0.8),
// Array(0.651, 0.682, 0.581, 0.741, 1, 0.798),
// Array(0.826, 0.75, 0.742, 0.8, 0.798, 1)):_*
// )