VBA Function for Eigen decomposition returns only positive eigenvalues

thx for feedback hagman, much appreciated...

would you have suggestion as to how I may modify/amend the code?

cheers

That would require me do dig deeper into your code.
I wonder a bit about lines 59/60:
I find
newA(i,j)^2 + newA(i,k)^2 = (A(i,j) * Cos_ + A(i,k) * Sin_)^2 + (-A(i,j) * Sin_ + A(i,k) * Cos_)^2
= A(i,j)^2 Cos^2 + 2 Sin Cos A(i,j) A(i,k) + A(i,k)^2 Sin^2  +  A(i,j)^2 Sin^2 - 2 Sin Cos A(i,) A(I,k) + A(i,k)^2 Cos^2
= A(i,j)^2 + A(i,k)^2
Hence if Application.SumSq(A) in line 67 does what its name suggests the test in line 68 should practically always evaluate to true, I think

Anyway, something like the following might help when inserted between lines 89 and 90:

If Ematrix(j,1) > 0 Then
 ' Find biggest component
 m = 2
 For i=3 to p+1
  If Abs(A(i,j)) > Abs(A(m,j)) Then 
    m = i
  End If
 Next i
 ' Calculate m'th component of M and jth eigenvector
 tmp = 0
 For i=1 to p
  tmp = tmp + M(m,i) * A(i,j)
 Next i
 ' recalculate eigenvalue as quotient output/input
 Ematrix(j,1) = tmp / A(m,j)
End If
                                              

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:

Thanks a lot for the follow-up Hageman...is there anything I can do to make you dig deeper into that code ;-)

I inserted your code btwn lines 89/90 (i.e. between the two next statements, pls see below); however, I know get an error msg 'subscript out of range')

Also, I noticed you set m= 2  & i = 3 to p+1....Wouldny p+1 cause the above error? Also, I should have mentioned that the input range maybe of variable size, hence in the smallest case, I may just have a 2x2 matrix....is there a more generic form of your fix which takes this into account?

Thanks a lot for your help - I really vm appreciate it!


Sub EIGEN_JKs()

Dim A() As Variant, Ematrix() As Double
Dim i As Long, j As Long, k As Long, iter As Long, p As Long
Dim den As Double, hold As Double, Sin_ As Double, num As Double
Dim Sin2 As Double, Cos2 As Double, Cos_ As Double, Test As Double
Dim Tan2 As Double, Cot2 As Double, tmp As Double
Const eps As Double = 1E-16
   
    On Error GoTo EndProc
   
    A = Range("L8:P12")
   
    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
       
        If Ematrix(j, 1) > 0 Then
 ' Find biggest component
 M = 2
 For i = 3 To p + 1
  If Abs(A(i, j)) > Abs(A(M, j)) Then
    M = i
  End If
 Next i
 ' Calculate m'th component of M and jth eigenvector
 tmp = 0
 For i = 1 To p
  tmp = tmp + M(M, i) * A(i, j)
 Next i
 ' recalculate eigenvalue as quotient output/input
 Ematrix(j, 1) = tmp / A(M, j)
End If
       
    Next j
       
       
    For i = 1 To p + 1 '
        For j = 1 To p
            Range("Z12").Offset(i, j) = Ematrix(j, i)
        Next j
    Next i
       
       
    'EIGEN_JK = Ematrix
   
    Exit Sub
EndProc:
    MsgBox prompt:="Error in function EIGEN_JK!" & vbCr & vbCr & _
        "Error: " & Err.Description & ".", Buttons:=48, _
        Title:="Run time error!"
End Sub