wilmott.com

HiI have below VBA code to determine the eigenvalues/vectors of a symmetric matrix. It seems that for some reason, the code only returns positive eigenvalues (although the absolute #s are correct)Eg for the following symmetric matrix1 22 1The code returns the Eigenvalues {3;1}, whereas the 'true' eigenvalues are {-1;3}The Eigenvectors though seem correct(0.707; 0.707) for +3(-0.707; 0.707) for 1Any idea how I may change the below code such that it gives me the correct sign?Thanks a lot!Public Function EIGEN_JK(ByRef M() As Variant) As VariantDim A() As Variant, Ematrix() As DoubleDim i As Long, j As Long, k As Long, iter As Long, p As LongDim den As Double, hold As Double, Sin_ As Double, num As DoubleDim Sin2 As Double, Cos2 As Double, Cos_ As Double, Test As DoubleDim Tan2 As Double, Cot2 As Double, tmp As DoubleConst eps As Double = 1E-16 On Error GoTo EndProc A = M p = UBound(A, 1) ReDim Ematrix(1 To p, 1 To p + 1) For iter = 1 To 15 'Orthogonalize pairs of columns in upper off diag For j = 1 To p - 1 For k = j + 1 To p den = 0# num = 0# 'Perform single plane rotation For i = 1 To p num = num + 2 * A(i, j) * A(i, k) den = den + (A(i, j) + A(i, k)) * _ (A(i, j) - A(i, k)) Next i 'Skip rotation if aij is zero and correct ordering If Abs(num) < eps And den >= 0 Then Exit For 'Perform Rotation If Abs(num) <= Abs(den) Then Tan2 = Abs(num) / Abs(den) Cos2 = 1 / Sqr(1 + Tan2 * Tan2) Sin2 = Tan2 * Cos2 Else Cot2 = Abs(den) / Abs(num) Sin2 = 1 / Sqr(1 + Cot2 * Cot2) Cos2 = Cot2 * Sin2 End If Cos_ = Sqr((1 + Cos2) / 2) Sin_ = Sin2 / (2 * Cos_) If den < 0 Then tmp = Cos_ Cos_ = Sin_ Sin_ = tmp End If Sin_ = Sgn(num) * Sin_ 'Rotate For i = 1 To p tmp = A(i, j) A(i, j) = tmp * Cos_ + A(i, k) * Sin_ A(i, k) = -tmp * Sin_ + A(i, k) * Cos_ Next i Next k Next j 'Test for convergence Test = Application.SumSq(A) If Abs(Test - hold) < eps And iter > 5 Then Exit For hold = Test Next iter If iter = 16 Then MsgBox "JK Iteration has not converged." 'Compute eigenvalues/eigenvectors For j = 1 To p 'Compute eigenvalues For k = 1 To p Ematrix(j, 1) = Ematrix(j, 1) + A(k, j) ^ 2 Next k Ematrix(j, 1) = Sqr(Ematrix(j, 1)) 'Normalize eigenvectors For i = 1 To p If Ematrix(j, 1) <= 0 Then Ematrix(i, j + 1) = 0 Else Ematrix(i, j + 1) = A(i, j) / Ematrix(j, 1) End If Next i Next j EIGEN_JK = Ematrix Exit Function EndProc: MsgBox prompt:="Error in function EIGEN_JK!" & vbCr & vbCr & _ "Error: " & Err.Description & ".", Buttons:=48, _ Title:="Run time error!"End Function

Check www.rnfc.org/ivey/ --> Math Section --> MATRIX_EIGEN --> There are different ways. Here it is the Jacobi version that give you the right values.Private Function PORT_FACTOR_EIGEN_VALUES_VECTORS_JACOBI_FUNC(ByRef DATA_RNG As Variant, _Optional ByVal OUTPUT As Integer = 0, _Optional ByVal tolerance As Double = 10 ^ -5)Dim i As LongDim j As LongDim ii As LongDim jj As LongDim kk As LongDim NSIZE As LongDim NROWS As LongDim NCOLUMNS As LongDim PI_VAL As DoubleDim MAX_VAL As DoubleDim RAD_VAL As DoubleDim TEMP_SUM As DoubleDim ANEXT_MATRIX As VariantDim VNEXT_MATRIX As VariantDim PTHIS_MATRIX As VariantDim ROTATION_ARR As VariantDim TEMP_MATRIX As VariantDim EIGENVECTORS_MATRIX As VariantDim EIGENVALUES_VECTOR As VariantDim DATA_MATRIX As Variant 'is rotated (using the P matrix)'until its off-diagonal elements are minimalOn Error GoTo ERROR_LABELPI_VAL = 3.14159265358979DATA_MATRIX = DATA_RNGNROWS = UBound(DATA_MATRIX, 1)NCOLUMNS = UBound(DATA_MATRIX, 2)NSIZE = Int(Sqr(NROWS * NCOLUMNS))If OUTPUT = 1 Then: GoTo 1983'------------------------------------------------------------------------------------------------'First Pass: Uses the Jacobi method to get the eigenvalues for a symmetric matrix'------------------------------------------------------------------------------------------------kk = 0TEMP_SUM = MATRIX_TRIANGULAR_UPPER_SUM_FUNC(DATA_MATRIX)Do While TEMP_SUM > tolerance GoSub P_ROTATION_LINE TEMP_SUM = MATRIX_TRIANGULAR_UPPER_SUM_FUNC(ANEXT_MATRIX) DATA_MATRIX = ANEXT_MATRIX kk = kk + 1LoopReDim EIGENVALUES_VECTOR(1 To NSIZE, 1 To 1)For i = 1 To NSIZE EIGENVALUES_VECTOR(i, 1) = DATA_MATRIX(i, i)Next iIf OUTPUT = 0 Then PORT_FACTOR_EIGEN_VALUES_VECTORS_JACOBI_FUNC = EIGENVALUES_VECTOR Exit FunctionEnd If1983:'------------------------------------------------------------------------------------------------'Second Pass: Uses the Jacobi method to get the eigenvectors for a symmetric matrix'Similar to eigenvalue function, but with additional V matrix updated with each rotation'------------------------------------------------------------------------------------------------kk = 0ReDim EIGENVECTORS_MATRIX(1 To NSIZE, 1 To NSIZE)For i = 1 To NSIZE: EIGENVECTORS_MATRIX(i, i) = 1: Next i 'Identity MatrixTEMP_SUM = MATRIX_TRIANGULAR_UPPER_SUM_FUNC(DATA_MATRIX)Do While TEMP_SUM > tolerance GoSub P_ROTATION_LINE GoSub V_ROTATION_LINE TEMP_SUM = MATRIX_TRIANGULAR_UPPER_SUM_FUNC(ANEXT_MATRIX) DATA_MATRIX = ANEXT_MATRIX EIGENVECTORS_MATRIX = VNEXT_MATRIX kk = kk + 1LoopEIGENVECTORS_MATRIX = VNEXT_MATRIXIf OUTPUT = 1 Then PORT_FACTOR_EIGEN_VALUES_VECTORS_JACOBI_FUNC = EIGENVECTORS_MATRIX Exit FunctionEnd If'------------------------------------------------------------------------------------------------PORT_FACTOR_EIGEN_VALUES_VECTORS_JACOBI_FUNC = Array(EIGENVALUES_VECTOR, EIGENVECTORS_MATRIX)'------------------------------------------------------------------------------------------------Exit Function'------------------------------------------------------------------------------------------------P_ROTATION_LINE:'Returns Anext matrix, updated using the P rotation matrix'------------------------------------------------------------------------------------------------ GoSub A_ROTATION_LINE GoSub R_ROTATION_LINE ANEXT_MATRIX = MMULT_FUNC(MATRIX_TRANSPOSE_FUNC(PTHIS_MATRIX), _ MMULT_FUNC(DATA_MATRIX, PTHIS_MATRIX, 70), 70)'------------------------------------------------------------------------------------------------Return'------------------------------------------------------------------------------------------------'------------------------------------------------------------------------------------------------V_ROTATION_LINE:'Returns Vnext matrix; keeps track of the eigenvectors during the rotations'------------------------------------------------------------------------------------------------ GoSub A_ROTATION_LINE GoSub R_ROTATION_LINE VNEXT_MATRIX = MMULT_FUNC(EIGENVECTORS_MATRIX, PTHIS_MATRIX, 70)'------------------------------------------------------------------------------------------------Return'------------------------------------------------------------------------------------------------'------------------------------------------------------------------------------------------------A_ROTATION_LINE: 'Returns vector containing the row and column vectors and'the angle of rotation for the P matrix'------------------------------------------------------------------------------------------------ ReDim TEMP_MATRIX(1 To NSIZE, 1 To NSIZE) MAX_VAL = -1 ii = -1 jj = -1 For i = 1 To NSIZE For j = i + 1 To NSIZE TEMP_MATRIX(i, j) = Abs(DATA_MATRIX(i, j)) If TEMP_MATRIX(i, j) > MAX_VAL Then MAX_VAL = TEMP_MATRIX(i, j) ii = i jj = j End If Next j Next i If DATA_MATRIX(ii, ii) = DATA_MATRIX(jj, jj) Then RAD_VAL = 0.25 * PI_VAL * Sgn(DATA_MATRIX(ii, jj)) Else RAD_VAL = 0.5 * Atn(2 * DATA_MATRIX(ii, jj) / (DATA_MATRIX(ii, ii) - DATA_MATRIX(jj, jj))) End If ROTATION_ARR = Array(ii, jj, RAD_VAL)'------------------------------------------------------------------------------------------------Return'------------------------------------------------------------------------------------------------'------------------------------------------------------------------------------------------------R_ROTATION_LINE:'Returns the rotation PTHIS_MATRIX matrix'------------------------------------------------------------------------------------------------ ReDim PTHIS_MATRIX(1 To NSIZE, 1 To NSIZE) For i = 1 To NSIZE: PTHIS_MATRIX(i, i) = 1: Next i 'Identity Matrix PTHIS_MATRIX(ROTATION_ARR(1), ROTATION_ARR(1)) = Cos(ROTATION_ARR(3)) PTHIS_MATRIX(ROTATION_ARR(2), ROTATION_ARR(1)) = Sin(ROTATION_ARR(3)) PTHIS_MATRIX(ROTATION_ARR(1), ROTATION_ARR(2)) = -Sin(ROTATION_ARR(3)) PTHIS_MATRIX(ROTATION_ARR(2), ROTATION_ARR(2)) = Cos(ROTATION_ARR(3))'------------------------------------------------------------------------------------------------Return'------------------------------------------------------------------------------------------------ERROR_LABEL:PORT_FACTOR_EIGEN_VALUES_VECTORS_JACOBI_FUNC = Err.numberEnd FunctionRafael Nicolas Fermin

Even better. Full matrix package for excel with source:Matrix 2.3 (http://digilander.libero.it/foxes/SoftwareDownload.htm)Cheers,Ram

rmeenaks,I totally agree. That's probably the best out there. I used it and continue to use it.Good to see that someone knows it.

The problem with the EIGEN_JK function is that it only works for positive definite matrices. If in the end of the function'Compute eigenvaluesFor k = 1 To pEmatrix(j, 1) = Ematrix(j, 1) + A(k, j) ^ 2Next kEmatrix(j, 1) = Sqr(Ematrix(j, 1))is replaced with'Compute eigenvaluesFor i = 1 To pFor j = 1 To pEmatrix(k, 1) = Ematrix(k, 1) + A(i, k) * A(j, k) * M(i, j)Next jNext iIf Ematrix(k, 1) > 0 ThenEmatrix(k, 1) = Ematrix(k, 1) ^ (1 / 3)ElseEmatrix(k, 1) = -(-Ematrix(k, 1)) ^ (1 / 3)End Ifthe function produces the correct eigenvalues also for non-definite matricesBr

I dont know why but the EIGEN_JK gives me an error.does the input neds to be positive definite? Moreover what should be input of this - a matrix range in excel with Ctr+Shft enter and output being another matrix area?

If e is an eigenvector of A with eigenvalue v thenA e = v ebut alsoA (-e) = v (-e)therefore the sign is arbitrary

Hiapologies for warming up this thread, but I was trying the apply the function introduced in the initial thread - Public Function EIGEN_JK(ByRef M() As Variant) As Variant - once more.In order to address the issue that the function fails to return negative eigenvalues in case of non psd matrices, I tried out Qazzt's suggestion, but it doesn't seem to work...Was wondering whether someone may have alternative suggestionscheersThe problem with the EIGEN_JK function is that it only works for positive definite matrices. If in the end of the function'Compute eigenvaluesFor k = 1 To pEmatrix(j, 1) = Ematrix(j, 1) + A(k, j) ^ 2Next kEmatrix(j, 1) = Sqr(Ematrix(j, 1))is replaced with'Compute eigenvaluesFor i = 1 To pFor j = 1 To pEmatrix(k, 1) = Ematrix(k, 1) + A(i, k) * A(j, k) * M(i, j)Next jNext iIf Ematrix(k, 1) > 0 ThenEmatrix(k, 1) = Ematrix(k, 1) ^ (1 / 3)ElseEmatrix(k, 1) = -(-Ematrix(k, 1)) ^ (1 / 3)End IfText

QuoteOriginally posted by: rmeenaksEven better. Full matrix package for excel with source:Matrix 2.3 (http://digilander.libero.it/foxes/SoftwareDownload.htm)Cheers,RamHi, gents!Seems that this link is dead, could anyone provide this fantastic excel package?

Well, it is documented what to do in page 272 of Kaiser's 1971 paper, The JK Method... this is where you'll find what Qazzt was trying to mention.But anyway it's correct, to replace the code as Qazzt said... just you need to do some reworking of the indices. The full (replacement) code is as follows:... above is as per original... 'Compute eigenvalues/eigenvectors For k = 1 To p For i = 1 To p For j = 1 To p Ematrix(k, 1) = Ematrix(k, 1) + A(i, k) * A(j, k) * m(i, j) Next j Next i If Ematrix(k, 1) > 0 Then Ematrix(k, 1) = Ematrix(k, 1) ^ (1 / 3) Else Ematrix(k, 1) = -(-Ematrix(k, 1)) ^ (1 / 3) End If 'Normalize eigenvectors For i = 1 To p Ematrix(i, k + 1) = A(i, k) / Ematrix(k, 1) Next i Next k EIGEN_JK2 = Ematrix Exit Function