jacobinomax.m

jacobi algorithm matlab code

%%%%%%%%%%%%%%%%%%%JACOBI ALGORITHM WITHOUT SORTING%%%%%%%%%%%%%%%%%%%%%%%

%%INPUT
%%randMatrix is the matrix we wish to diagonalize. It's only random for the
% first interation
%%k is our current iteration
%%A is the matrix previous to randMatrix 

%%OUTPUT
%%D is the more diagonal matrix
%%bk is the off(D)
%%k is the number of iterations

function [D, bk, k, originalMatrix] = jacobiNoMax(randMatrix, k, A)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%Creating the random symmetric matrix [START]%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if (size(randMatrix)==[0,0])
        
        %randMatrix is our random 5 x 5 symmetric matrix and is initialized
        %to be completely filled with 0's
        randMatrix = zeros(5,5);
        
        %If the row is equal to the column then we are on the diagonal
        %and thus the values needn't be mirrored. If we are not on
        %the diagonal, then we need to mirror the random value on the
        %other side of the diagonal to create a symmetric matrix. When 
        %a random value is stored in randMatrix(row, column), it is also 
        %stored in (column, row), thus "mirroring" the value. 
        %
        %We will loop through the x portion of the randMatrix:
        %      | x 0 0 0 0 |
        %      | x x 0 0 0 |
        %      | x x x 0 0 |
        %      | x x x x 0 |
        %      | x x x x x |
        %and mirror the off-diagonal elements. The diagonal elements will a
        %sense be "mirrored" too, but mirrored with its own position. 
        %Ex: row=2, column=2, randNumber is some random value
        %    then randMatrix (row, col) = randNumber;
        %         randMatrix (col, row) = randNumber;
        %    is   randMatrix (2,2) = randNumber;
        %         randMatrix (2,2) = randNumber;
        
        for (row = 1:5)
            for (col = 1:row)
                  %randNumber is a random number on the interval [0,10]  
                  randNumber = floor(rand*10);
                  %change the value at randMatrix(row, col) to the random number
                  randMatrix (row, col) = randNumber;
                  %mirror the value at randMatrix(col, row)
                  randMatrix (col, row) = randNumber;
            end
        end
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%Creating the random symmetric matrix randMatrix [END]%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

originalMatrix = randMatrix;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%Randomly selecting an off-diagonal entry and the 2 x 2 matrix%%%%%%%% 
%%%%%%containing that entry [START]%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %B will our 2 x 2 matrix
    B = zeros(2,2);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%ONLY PORTION DIFFERENT FROM JACOBI ALGORITHM WITH SORTING%%%%%%%%%    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
  
    %Instead of looking for the off-diagonal entry with the largest
    %magnitude. We instead randomly choose an off diagonal entry.
    maxRow = 0;
    maxCol = 0;
    
    %Since our modified rand() will generate random numbers from 0 to 5,
    %we need to loop until neither numbers are 0 because matlab begins
    %indexing at 1 instead of 0 like most other languages and until the
    %two numbers are not equal, because we wish for an off-diagonal entry
    %rather than a diagonal entry.
    %%
    while (maxRow == 0)||(maxCol == 0)||(maxRow == maxCol)
        maxRow = abs(floor(rand()*10)-4);
        maxCol = abs(floor(rand()*10)-4);
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %So we've found the entry. Now we need to use it to find a 2 x 2
    %symmetric matrix that we can diagonalize. The matrix basically follows
    %this formula which should hold true any n x n symmetric matrix where 
    %n >=3:
    %B = | randMatrix(maxRow, maxRow)  randMatrix(maxRow, maxCol) |
    %    | randMatrix(maxRow, maxCol)  randMatrix(maxCol, maxCol) | 
    B(1,2) = randMatrix(maxRow, maxCol);
    B(2,1) = randMatrix(maxRow, maxCol);
    B(1,1) = randMatrix(maxRow, maxRow);
    B(2,2) = randMatrix(maxCol, maxCol);
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%Randomly selecting an off-diagonal entry and the 2 x 2 matrix%%%%%%%% 
%%%%%%containing that entry [END]%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%Diagonalizing 2 x 2 symmetric matrix B [START]%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % a, b, d will be used the calculation of the eigenvalues where:
    % B = |a b |
    %     |b d |
    a = B(1,1);
    b = B(1,2);
    d = B(2,2);
    I = [1 0; 0 1];
    
    %In order to diagonalize the matrix, we set the det(B-(mu*I)) equal to 
    %0 and solve for mu:
    %    0 = det | a - mu  b    |
    %            | b       d-mu |
    %   0 = (a-mu)*(d-mu) - b^2
    %   0 = mu^2 - (d+a)mu + (ad-b^2)
    %   In this form, we can use the quadatric formula to solve for mu:
    %       mu+ = myPlus  = ((d+a) + sqrt((d-a)^2 + 4*b^2))/2
    %       mu- = muMinus =((d+a) - sqrt((d-a)^2 + 4*b^2))/2
    %But we also need to check if b is equal to 0, because if it is then
    %this matrix is already diagonalized. We can eliminate needless
    %calculations and avoid incidents that will cause our matrix to be 
    %filled with NaN (not a number) values. ie divide by zero or squareroot
    %of a negative number.
    if (b ~= 0)
        muPlus = ((d+a) + sqrt((d-a).^2 + 4.*b.^2))./2;
        muMinus = ((d+a) - sqrt((d-a).^2 + 4.*b.^2))./2;
        %Now we plug in the eigenvalues where u1 and u2 are the results of that
        %plugging in, though it's not really necessary that we calculate the
        %u2
        u1 = B - (muPlus*I);
        if (u1 == NaN)
            u1 = B - (muMinus*I);
        end
        u2 = B - (muMinus*I);

        %u1R1 is the first column of our soon to be normalized matrix u. u1R1
        %is obtained from the first row of the u1 with the appropriate sign
        %flipping. Ex:  u1 = |a b|
        %                    |c d|
        %               u1R1 = |-b|
        %                      | a|
        %but technically we could skip the skip of actually looking for the
        %eigenvalues is the sum of the first row is equal to the sum of the
        %second. Ex:    | 2 1 | = 3
        %               | 1 2 | = 3
        %We can assume that |1| to be our u1R1.
        %                   |1|
        u1R1 = [-u1(1,2) u1(1,1)]';
        if (sum(B(1,:))==sum(B(2,:)))
            u1R1 = [1 1]';
        end

        %u2R1 is the second column of matrix u and is obtained from u1R1 again
        %with the appropriate flipping. 
        %           Ex: u1R1 = |a|
        %                      |b|
        %               u1R2 = |-b|
        %                      | a|
        u2R1 = [-u1R1(2) u1R1(1)]';

        %u1Length is the 1/||u1R1||. The length for u2R1 is the same as u1R1 so
        %there is no need to calculate it again.
        u1Length = 1/(sqrt(u1R1(1)^2 + u1R1(2)^2));

        %The first column of matrix u is u1R1 and the second is u2R1. u is
        %normalized because its components are both normalized.
         u = u1Length*[u1R1 u2R1];


        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%Creating the Givens Matrix [START]%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%G will be our Givens matrix, first filled with 0's 
            G = zeros(5,5);

            %We place entries from the u into G that correspond to the entries
            %that were taken from the original matrix randMatrix
            G(maxRow, maxCol) = u(1,2);
            G(maxCol, maxRow) = u(2,1);
            G(maxRow, maxRow) = u(1,1);
            G(maxCol, maxCol) = u(2,2);

            %We loop through G placing 1's on the diagonal where the value
            %is 0.
            for row=1:5
                for col=1:5
                    if (row==col)
                        if (G(row,col)==0)
                            G(row,col)=1;
                        end
                    end
                end
            end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%Creating the Givens Matrix [END]%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        %D is our final more-diagonal matrix that we return
        D = G'*randMatrix*G;
    else
        D = randMatrix;
    end
    %bk is the value of the OFF(D), where OFF(x) is the sum of the square of 
    %each off-diagonal elements of x 
    bk = 0;
    for (row = 2:5)
            for (col = 1:row-1)
                bk = bk + 2*(D(row, col))^2;
            end
    end

    %increment k
    k=k+1;
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%Diagonalizing 2 x 2 symmetric matrix B [END]%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%