Ternary Search Tree in C/C++

Task is to store an array of integers in a ternary search tree. Each element of the array can take one of the three values: 0, 1 or 2.

The C/C++ code is available on my GITLAB page.  The resulting tree could be plotted using Graphviz software using a dot file.  This code also contains code for generating a dot file. Once the dot file is obtained. Use the following command to generate the tree as shown below:

$ dot -Tps tst_data.dot -o tst_data.ps
$ gv tst_data.ps

Data Generation for Tic-Tac-Toe Game : Example of Recursion in C++

Problem Statement: I need to generate data set for a 3x3 Tic-Tac-Toe Game. In this game, there are 9 Cells and each can take one of the values in the set {0, 1, 2}. Where 0 denotes and empty cell while 1 and 2 correspond to each of the two players.  So, we need to generate a 9x1 state vector where each element can have one of three possible states. Total number of states is 3^9 = 19683.


I have written a generic program that can generate such datasets irrespective of their dimension or possible state values:

We use a function to remove redundant vectors.

#include <iostream>
#include <fstream>
#include <cmath>

using namespace std;

#define N 3
#define M 3
#define MAX 20000
int Data[MAX][N];

std::ofstream f1("dataset2.txt");
//-------------------------------------------

bool find_match(int state[N], int rowcnt)
{
  for(int q = 0; q < rowcnt; q++)
  {
    int simcnt = 0;
    for(int r = 0; r < N; r++)
    {
      if(Data[q][r] == state[r])
        simcnt++;
    }

    if(simcnt == N)
      return true;
  }

  return false;
}

//------------------------------------

void add_row(int state[N], int rowcnt)
{

  for(int r = 0; r < N; r++)
    Data[rowcnt][r] = state[r];
}

//-----------------------------------------

void genstate(int state[], int k, int &rowcnt)//, char *filename)
{

  bool mFlag = false;
  for(int l = k; l < N; l++)
  {
    for(int j = 0; j < M; j++)
    {
      state[k] = j;

      mFlag = find_match(state, rowcnt);
      if(mFlag == false)
      {
        add_row(state, rowcnt);
        rowcnt++;
        for(int p = 0; p < N; p++)
          f1 << state[p] << "\t";
        f1 << endl;
      }
      genstate(state, k+1,rowcnt);

    }           
  }

 // f1.close();
}

//-------------------------

int main()
{
  int state[N];
  for(int i = 0; i < N; i++)
    state[i] = 0;

  int rowcnt = 0;

  genstate(state, 0, rowcnt);//, "dataset.txt");

  f1.close();

  return 0;
}


       
 

Systematic Resampling

The idea is simple, but I don't know why I took so much of time to understand it. I am still not sure if I got it right. Anyway, let me explain what I understood.

Please refer to this page to understand the basics of systematic sampling. The starting point is chosen at random and then the other points are selected at regular intervals. 

Let us assume that we have population X[N] with N elements. The population has an associated weight vector W[N]. We need to sample Y[n] of n elements out of X using systematic sampling. The algorithm that I could understand from Section 8.6 of [1] is as follows:

1. Compute the sum of all weights $ \displaystyle S = \sum_i^N W_i $
2. Compute the selection interval I = S/n ; 
3. To initiate the random selection process, select a uniform random number R in the open interval (0, I].
4. The n samples are given as : R, R+I, R+2I, ... , R+(n-1)I

void systematic_resampling1(gsl_rng *r, double *dest, size_t m, double *src, size_t n, double *w)

{

  if(m > n)

  {

    cerr << "choose m <= n" << endl;

    exit(-1);

  }

  float sum_w = 0.0;

  for(size_t i = 0; i < n; ++i)

    sum_w += w[i];

  // selection interval

  float I = sum_w / m ;  

  //generate a random number between 0 and I

  float u = gsl_rng_uniform_pos(r) * I;

  size_t j = 0;

  while(j < m)

  {            

    float C = 0.0;

    for(size_t i = 0; i < n; ++i)

    {

      C = C + w[i];      

      if(u <= C)

      {

        dest[j] = src[i];

        j = j + 1;

        u = u + I;

      }

    }

  }

}

The other Algorithm [2], that seems to achieve the same is as follows:

void systematic_resampling2(gsl_rng *r, double *dest, size_t m, double *src, size_t n, double *w)
{
  if(m > n)
  {
    exit(-1);
  }
  float sum_w = 0.0;
  for(size_t i = 0; i < n; ++i)
    sum_w += w[i];

  //generate a random number between 0 and 1
  float u = gsl_rng_uniform_pos(r);

  int j = 0;
  float C = 0.0;
  for(size_t i = 0; i < n; ++i)
  {
    C = C + w[i]/sum_w;      
    if(C > 1.0)
      C = 1.0;

    while((u+j)/m <= C)
    {
      dest[j] = src[i];
      j = j + 1;
    }
  }
}

Reference:

[1]  Kirk M. Wolter, "Introduction to Variance Estimation", Statistics for Social and Behavioral Sciences series, Springer, 2nd Edition, 2007.

[2] Murray Pollock, "Introduction to Particle Filtering", Algorithms & Computationally Intensive Inference Reading Group Discussion Notes, 2010. 

Weighted Random Sampling With Replacement in C

The problem is to select k items with-replacement from a list of n items with a probability that depends on the weights associated with the items in the original list.

Input:  x[n], w[n]
Output: x[n]   (x is overwritten)

void wrswr(gsl_rng *r, std::vector &x, const std::vector &w)
{
  if(x.size() != w.size())
  {
    cerr << "w and x must have same length ..." << endl;
    exit(-1);
  }
  else
  {
    size_t N = (size_t)x.size();
    std::vector tempx(N);

    double sum_w = 0.0;
    for(size_t i = 0; i < N; ++i)
    {
      sum_w += w.at(i);
      tempx.at(i) = x.at(i);  // copy of source xk
    }

    for(size_t i = 0; i < N; ++i)
    {
      double u = gsl_rng_uniform(r);

      size_t j = 0;
      double c = w[0] / sum_w;
      while (u > c && j < N)
      {
        j = j + 1;
        c = c + w[j] / sum_w;
      }
      x.at(i) = tempx.at(j); 
    }
  }
}         

Weighted Random Sampling

You need to compile it using GSL library. I implemented the idea available on this post.  For a normal random sampling with replacement algorithm without weights is available in GSL itself. Refer to this page for details on "Shuffling and Sampling" with GSL. 

Latex on Blogger

I have been on blogger for a long time now. I don't want to move over to Wordpress which is known to provide better support for mathematical equations. However, there are many like me and people have found solutions. There are many ways to do it. But  this solution worked for me.

Just a test :

\int_{0}^{1}\frac{x^{4}\left(1-x\right)^{4}}{1+x^{2}}dx
=\frac{22}{7}-\pi

Another solution is available here. It requires you to install a script into your browser (chrome or firefox). After that it shows two extra menu buttons on your blogger editor : Latex and UnLatex. A latex code put in between a pair of  'double dollars' converted into an image equation. You can undo the operation by clicking on UnLatex. Another test:

Bias Vs Variance Dilemma

Bias versus Variance dilemma is almost universal in all kinds of data modelling methods. As per Wikipedia, variance of a random variable is the deviation squared of that variable from its expected value. In other words, it is a measure of variation within the values of the variable across all possible values along with their probabilities. The bias of an estimator is the difference between the estimator's expected value and the true value of the parameter being estimated. A very good article on this topic is available here. Without repeating much of the content, I would simply highlight the key points which would make things easier in understanding this topic.

1. Var(x) = E(x^2) - [E(x)]^2, where E(.) is the expectation value.  This can be rewritten as
   
   
   E(x^2) = Var(x) + [E(x)]^2
   If we replace x by e (approximation error of an estimator), we can rewrite above equation as
   
   E(e^2) = Var(e) + [E(e)]^2
   MSE = Var(e) + Bias^2
   
   Hence we can see that for a desired mean square error, there is trade-off between the variance and the bias. If one increases, the other decreases.
   
   
2. The complexity of an estimator model is expressed in terms of the number of parameters. The effect of complexity on the bias and the variance of the estimator is explained here. In brief, it can be said that
1. Low Complexity leads to low variance but large bias.
2. Highly complex model leads to low bias but large variance.
   
   
3. Large variance implies that the estimator is too sensitive to the data set. Hence the model has a low generalization capability. Excellent performance on design (or training) data, but poor performance on test data. This is the case of overfitting. Large variance is observed in case of complex models with large number of parameters (over-parameterization). Note that because of high complexity, the model fits a  the design data set very well. Hence, the error at individual points (in the data set) is very low and hence the model has a low bias.
   
   
4. Large bias implies that the model is too simple and hence very few data points lie on the regression curve. Hence, the error at individual points (in a given data set) are high leading to large a MSE. But since very few points participate in the model formation, the performance does not differ on different data tests. Hence, it has a low variance. The performance remains same over design and test data sets.

Removing some rows from a matrix, Matlab

>> A = magic(4)
A =
    16     2     3    13
     5    11    10     8
     9     7     6    12
     4    14    15     1
>> b = [ 2     3];
>> A(b,:) = []
A =
    16     2     3    13
     4    14    15     1

More on manipulating 1-D array

Matlab Path

This is how I add a folder to the matlab path on command window:



if(~exist('varycolor'))
    addpath(genpath('~/MatlabAddon'));
end


The if condition checks for the existence of the function 'varycolor'. If it does not exists, it loads the folder where this file is present.

Deleting empty rows in a matrix

>> a = [1,2,3;4,5,6;0,7,8;0,0,0;0,0,0]

a =
     1     2     3
     4     5     6
     0     7     8
     0     0     0
     0     0     0

>> b = a(any(a,2),:)

b =

     1     2     3
     4     5     6
     0     7     8
>> b = a(all(a,2),:)

b =

     1     2     3
     4     5     6

K-Medoid Algorithm in Matlab

According to Wikipedia, k-medoids algorithm is a clustering algorithm related to the k-means algorithm. In contrast to k-means algorithm, k-medoids chooses data points as centres.  I wrote a Matlab program for implementing this algorithm. I am posting it here in the hope that people may find it useful.

function [IDX, Cluster, Err] = kmedoid2(data, NC, maxIter, varargin)
% Performs clustering using Partition Around Medoids (PAM) algorithm
% Initial clusters must be selected from the available data points

% Usage
% [IDX,C,cost] = kmedoid(data, NC, maxIter, [init_cluster])
% 
% Input : 
%       data - input data matrix
%       NC - number of clusters
%       maxIter - Maximum number of iterations
%       init_cluster - Initial cluster centers (optional)
%       init_index - Index of initial clusters
%
% Output :
%       IDX - data point index which became cluster centres
%       C -  Cluster centers
%       Err - cost associated with the clusters for all iterations
% ------------------------------------------

% Size of data
dsize = size(data);

% Number of data points
L = dsize(1);

% No. of features
NF = dsize(2); 

% Cluster size
csize = [NC, NF];

if (L < NC)
    error('Too few data points for making k clusters');
end

if(nargin > 5)
    error('Usage : [IDX,C,error] = kmedoid(data, NC, maxIter, [init_cluster], [init_idx])');
elseif(nargin == 5)
    vsize1 = size(varargin{1});
    vsize2 = length(varargin{2});
    
    if(isequal(vsize1,csize))
        Cluster = varargin{1};
    else
        error('Incorrect size for initial cluster');
    end
    
    if(vsize2 == NC)
        IDX = varargin{2};
    else
        error('Incorrect size for initial cluster index');
    end
elseif(nargin == 4)
   error('You must provide two optional arguments: init_cluster, init_idx');
elseif(nargin == 3) % no initial cluster provided
    IDX = randint(NC,1,L)+1;
    Cluster = data(IDX,:); % Initialize the initial clusters randomly
else
    display('Usage : [IDX,C] = kmedoid(data, NC, maxIter, [init_cluster], [init_idx]');
    error('Function takes at least 3 arguments');
end

%Array for storing cost for each iteration
Err = zeros(maxIter,1);

Z = zeros(NC, NF, L);
for i = 1:L
    Z(:,:,i) = repmat(data(i,:),NC,1);
end


FIRST = 1;
total_cost = 0;


for Iteration = 1:maxIter
    %disp(Iteration);
    
    if(FIRST) 
        % First time, compute the cost associated with the initial cluster
        
        C = repmat(Cluster,[1,1,L]);
        B = sqrt(sum((abs(Z-C).^2),2)); %Euclidean distance
        B1 = squeeze(B);
        [Bmin,Bidx] = min(B1,[],1);

        cost = zeros(1,NC);
        for k = 1:NC
            cost(k) = sum(Bmin(Bidx==k));
        end
        total_cost = sum(cost);
        
        FIRST = 0; % Reset the FIRST flag
        
    else % Not first time
        
        % change one cluster center randomly and see if the cost is
        % decreased. If yes, accept the change, otherwise reject the change

        while(1)
            Tidx = randint(1,1,L)+1; % find a random number betn 1 to L
            if(isempty(find(IDX==Tidx))) % Avoid previously selected centers
                break;
            end
        end
        % randomly change one cluster
        pos = randint(1,1,NC) + 1;
        OLD_IDX = IDX; % Preserve the old index list
        IDX(pos) = Tidx;

        %Assign cluster centers
        PCluster = Cluster;
        Cluster = data(IDX,:);

        % For each data point, find the cluster which is closest to it
        % and compute the cost obtained from this association

        C = repmat(Cluster,[1,1,L]);

        %B = sum(abs(Z-C),2); %Manhattan distance
        B = sqrt(sum((abs(Z-C).^2),2)); %Euclidean distance

        % Row - cluster (NC)
        % Column - each data point 1, 2, ... L
        B1 = squeeze(B);

        % For each data point, find the nearest cluster
        % size(Bmin) = 1xL
        [Bmin,Bidx] = min(B1,[],1);

        cost = zeros(1,NC);
        for k = 1:NC
            cost(k) = sum(Bmin(Bidx==k));
        end

        previous_cost = total_cost;
        total_cost = sum(cost);


        if(total_cost > previous_cost) % cost increases
            % Bad choice, restore the previous index list
            IDX = OLD_IDX;
            total_cost = previous_cost;
            Cluster = PCluster;
        end
    end
    Err(Iteration) = total_cost;
end

if(nargout == 0)
    disp(IDX);
elseif(nargout > 3)
    error('Function gives only 3 outputs');
end   
    
return  

Save this file as "kmedoid.m". Now we can create clusters for a given data set as follows :

[Idx, C,err] = kmedoid(data, 5, 1000);

Above command iterates for 1000 times to create 5 clusters for a given data.

Idx - the indices of data points selected as cluster centres.
C - cluster centres selected from the data sets
err - cost for each iteration. It can plotted to see if the error is decreasing.

Cycle through Markers in Matlab Plot

figure(2);
markers = 'ox+*sdv^<>ph';
colorset = varycolor(NC);

hold all;
for i = 1 : NC
    plot(C(:,1), C(:,2), 'LineStyle', 'None', 'Marker', markers(i), 'Markersize', 10);
end

Avoid loops in Matlab

Following examples demonstrate how one can avoid for loops in MATLAB to speed up execution :

1. Vector assignment
   
   x = rand(1,100);
   %In spite of :
   for k=1:100
   y(k) = sin(x(k));
   end
   % We can use :
   y(:) = sin(x(:));
   
   
   
2. Subtracting a vector from each row/column of a matrix
   
   a = [1 2 3]
   b = magic(3);
   c = repmat(a,3,1) - b;
   
   
3. Use logical indexing to operate on a part of a matrix. Use "find" command wherever possible.
   
   a = 200:-10:50;
   % replace the value 120 by 0
   a(find(a == 120)) = 0;
   a =   200   190   180   170   160   150   140   130     0   110   100    90    80    70    60    50
   
   Another example :
   
   % [idx1, idx2] = find(angle < 0);
   % idx = [idx2 idx2];
   % for j = 1:length(idx)
   %     angle(idx(j,1),idx(j,2))=360+angle(idx(j,1),idx(j,2));
   % end
   can be replaced by following code :
   
   angle(angle<0) = angle(angle<0)+360;
   
   
   
4. Useful commands for matrices and vectors :
   
   sum, norm, diff, abs, max, min, sort
   
   
   

   Useful link in this direction

   http://www.ee.columbia.edu/~marios/matlab/matlab_tricks.html  

   
   

Index of an array element closest to an arbitrary value

>> x = 1:100:1000;

The index of the array which is closest to say 130 is given by

>> [v,idx] = min(abs(x-130));

The answer is

v = 29, idx = 2

Effect of Controllability and Observability On Stability

In this post, I would post the email conversation that I had with my friends over the topic "the effect of controllability and observability on stability". I am thankful to the student Asmini who actually raised this question.

My response :

In simple terms, controllability is the virtue of a system by which one can control the state of the system. In technical terms, controllability is the ability to be able to control each state variable or to be able to take it from one point in state space to another point. There is subtle difference between the terms "Controllability" and "Reachability". Don't worry about this second term at present. Both are more or less same. If a system is controllable, that means we can design a controller to control all states. That's good for us, is n't it?

On the other hand, observability is the virtue of system which enables us to observe all the states. I mean if a system is observable, we can use sensors to record all state variables. Note that this is not always the case. We may not have access to all states of a system. Again, it is good to have a completely observable system.

Now coming to your question, what is the effect of controllability and observability on stability?

Well, I think if a system is completely controllable or observable, then it is possible to design controllers which can stabilise the system. But if the system is either uncontrollable or unobservable, then it would be difficult to stabilise it. In some case, it would be impossible to design a controller to stabilize the system.

Indrani's Response:

I assume that we know the meaning of "stability of a system". If a system
dynamics is not stable we can make it stable by application of a proper
control input. On the other hand, controllability of a system implies
that it is possible to force the system to a particular state by
application of a control input. If a state is uncontrollable then no
input will be able to control that state. Thus if a system state is not
controllable as well as  the corresponding dynamics is not stable then it
will not be possible to make that state stable by applying a control
input.


Similarly if a state is not observable then the controller will not be
able to determine its behavior from the system output and hence not be
able to use that state to stabilize the system. Thus controllability and
observability of a system are two important properties which are to be
considered before designing a controller.

Prem's response:

Yes... Indrani's answer is correct... There is no correlation between stability and controllability. A controllable system can be either unstable or stable.

But controllability means that if the system is unstable then you can stabilize it using proper state feedback. Observability indicates whether we can deduct the states from the observed output so that states will be available for stabilization.

If some of the nodes are not observable and those nodes are already stable then also we can stabilize.

so we can classify into two cases.

i) An observable and controllable system can be stabilizable with state feedbacm.

ii) If some of the states are unobservable and those nodes are stable, then again we can design state feedback to design the stabilizing the controller with observable states.

I am still waiting for Awhan's response and I would update this post once I hear from him. But following thing can be concluded about this topic.

• Controllability and Observability properties are desirable for controller design. This is the first thing to be checked before designing a controller.
  
  
• There is no correlation between Controllability and Observability with Stability of the system. Stability is the property of the system (Plant + Controller). On the other hand controllability and Observability are the properties of Plant alone.
  

Matlab Plotting within a loop with different colors

Matlab by default has 7 colors. Following code demonstrates how one can invoke plot command within a for loop so that each time it plots a given data set in a different color. Don't worry about the variables that are used in this piece of code. You should notice that we can pass randomized color information to the plot command inside a for loop.


hold on;
for k = 1:length(idx)
ridx = find(idx == k);
color = rand(1,3);
for j = 1:length(ridx);
plot(logdata(:,1,ridx(j)),logdata(:,2,ridx(j)), ...
'Color',color,'LineWidth',2);
end
end
hold off;

Useful Octave Links

Some good tutorials for beginners are available at following links:

1. Tutorial 1
   
   
2. Tutorial 2
   
   
3. Octave 3.0 documentation
   
   
4. Now a GUI front end is also available for Octave. It is called Qt-Octave. Click here to see a demo.
   

Octave and C++

I have been looking for a vector matrix library which can be used in C/C++ directly. Don't you think it would be cool if you could write something like C = A * B where C, A and B are matrices. This is not like I was n't aware of existence of such libraries. For instance, sometime back I used a library called CVMLIB. This library provides very good interface for all kinds of matrices and vectors along with a number of linear algebra routines. I don't know why I gave up using that. More recently I came across a library called liboctave. This provides a similar interface with the promise that it is possible to call several octave routines directly from a C++ program. I found it simpler to use as compared to CVMLIB. But I must admit that CVMLIB is more complete as compared to liboctave. If I can learn using octave functions, then it would be really great. Presently I am using a mixture of gsl, lapack, blas routines to get my work done.

Least Square Solution

I always get confused with this topic. Awhan explained it to me very nicely. I would discuss this issue briefly so that I can refer back to it when ever needed.

A set of linear equations is given by

A x = y .................... (1)

where A is m x n dimensional constant matrix. 'x' is a n-dimensional vector and y is a m-dimensional output or measurement vector. The task is to find a solution 'x' that satisfies equation (1).


Solution to equation (1) exists if and only if y lies in the range space of A. In other words,

rank[ A y] = rank[A] = r ............ (2)


In case this rank condition is satisfied, we know that there exists at least one solution. General solution to equation (1) is given by

x = xp + \lambda * xh .............. (3)

xp is the particular solution and xh is the homogenous solution. Multiplying both sides of equation (2) by A, we get

Ax = A * xp + \lambda A * xh = A * xp = y

Since xh is a vector that lies in the null space of A, we have A * xh = 0. Dimension of null space is given by \gamma = n - rank(A). Whenever \gamma > 0, or we have a non-trivial null space, we will have multiple solutions.


Lets assume that we have a solution for the set of equations given by (1). Consider following cases :

• If m is less than n . Less number of equations and more number of variables. This is also known as under-determined set of equations. If A is of full rank, that is , rank(A) = m, there exists a null space of dimension (n-m) and hence there would be many solutions to equation (1). You have many solutions even when rank(A) = r .
  
• If m > n. More number of equations and less number of variables. This is also known as overdetermined set of equations. If A is of full rank, that is, rank(A) = n, then the dimension of null space = n - n = 0. Hence there exists a unique solution to the set of equations given by (1). If rank(A) = r < n , then the set of equations (1) would have infinitely many solutions as the null space becomes non-trivial.

Date: 21 September 2008 Sunday

For m >= n (overdetermined system), the least square solution is obtained by minimizing the following cost function:

V = || Ax-y || ^2

which gives x = inv(A'*A) * A' * y = pinv(A)*y. This is also known as moore-penrose pseudo-inverse solution.

One can use matlab's backslash operator to compute the least square solution.

x = A\y

If A is not full rank, then A\b will generate an error message, and then
a least-squares solution will be returned.

-----------------

For under-determined system (m is less than or equal to n), the minimum-norm is solution is given by

x = A' * inv(A*A') * y

You can also use the pseudo-inverse function pinv(), which computes the pseudo-inverse,

x = pinv(A) * y;

when A is full rank and fat or square. (The pseudo-inverse is also defined when A is not full
rank, but it's not given by the formula above.)


Warning: If A is fat and full rank, then A\y gives a solution of the Ax = y, not necessarily the
least-norm solution.

The minimum-norm solution for least square could be derived by using Lagrange multiplier method:

Consider the problem of  minimizing ||x||^2 subject to Ax = y. For this, the Lagrange function may be written as

L(x, \lambda) = 0.5 x^x + \lambda^T * (Ax-y)

dL/dx = x^T + \lambda^T * A = 0.

This gives, x = A^T = \lambda and \lambda = inv (A*A^T) * y. Substituting the \lambda into the former, we get

x = A^T * inv(A*A^T) * y = pinv(A) * y

Non-minimum phase systems and Bode plot

Consider the bode plot for the system shown below:

What can be said about this system:

• The magnitude plot has a slope of -40 dB/decade. Its predominantly a second order system with double pole at origin.
  
• It has a non-minimum phase characteristic with negative phase margin.

Now can we say that the system is BIBO unstable. One of my undergrad student pointed out this to me. He said that its an unstable system because it seemed to have a double pole at origin. Infact another student approximated this bode plot with following transfer function:

s = tf('s');
g = 507*(s+7)/(s^2*(s^2+601*s+600));
bode(g);


However, we need to keep following points in mind:

• Non-minimum phase systems are not always unstable system.
• Negative phase margin does not always imply instability.
• For non-minimum phase systems with stable poles, we can't apply the approximations that hold for a standard second order system.
  

Controllability and Condition Number

Consider following System (A, B) given as follows:

A =
-205.5237 198.3209 7.2028 0 -0.3256 0 0
198.3209 -205.5237 0 7.2028 0 0 0
0.0646 0 -0.0646 0 0 0 0
0 0.0646 0 -0.0646 0 0 0
0 0 0 0 0 0 0
463.7826 0 0 0 0 -43.4783 0
0 463.7826 0 0 0 0 -43.4783

B =

0
0
0
0
0.5600
0
0


Lets check the controllability of this system:

>> rank(ctrb(A,B)) = 4

You get a different result when you apply PBH (Popov-Belevitch-Hautus) test. In this test we check for the rank of [\lambda*I-A, B] matrix for each eigenvalue lambda. Following code (written by Gopal) performs the PBH test


% %%Controllabilaty test
Eigen_A=eig(A);
Data_Con = zeros(7,2);
for i=1:length(Eigen_A)
S=Eigen_A(i)*eye(length(A) );
D=S-A;
Control_A=[D B];%Formation of controllability matrix.
Rank_Con = rank(Control_A);
Data_Con(i,:)= [Eigen_A(i), Rank_Con];
end
Data_Con
%%%%%%%%%%%%%%%%%%
This gives

Data_Con =

-43.4783 6.0000
-43.4783 6.0000
-403.8457 7.0000
-7.2674 7.0000
0.0000 7.0000
-0.0634 7.0000
0 7.0000

where first column gives the eigen values and the second gives the corresponding rank.

As far as I know the controllability information obtained in either way should be equivalent. That means both methods should give same information about the controllability.

On little analysis I found that the matrix A is highly ill-conditioned with a condition number

>> cond(A) = 4.3373e+20

Could this be the reason for this anomoly. In other words, you can not rely on the controllability tests as it is ill-conditioned. Gopal tells me that this system pertains to a real system (related to a nuclear reactor) . In that case, how do you deal with such a system. Is there any process by which we can make it well-conditioned so that one can carry out computer simulations with good accuracy.