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.

Stiff ODEs

Methods for solving differential equations may be divided into two classes namely explicit and implicit methods. Explicit methods calculate the state of a system at a later time from the state of the system at the current time, while an implicit method finds it by solving an equation involving both the current state of the system and the later one. Implicit methods require an extra computation and they can be much harder to implement. More information is available here. Implicit methods are used because many problems arising in real life are stiff for which the use of an explicit method requires impractically small time steps Δt to keep the error in the result bounded.

According to wikipedia, a stiff equation is a differential equation for which certain numerical methods for solving the equation are numerically unstable unless the step size is taken to be extremely small. The main idea is that the equation includes some terms that can lead to rapid variation in the solution.
We have seen that for solving ODE numerically the equation must be stable. This for example implies that the Jacobian of the RHS of the equation \dot{x} = f(x,t) must be a stable matrix. But with the Jacobian can be associated the Lipschitz constant, related with the norm of the Jacobian. Roughly a stiff problem is the combination of stability, with big eigenvalues (negative real part) of the Jacobian , implying a big norm, hence a big Lipschitz constant. This imply that the process described by the ODE contains components operating on different time scale.

For stiff problems, certain implicit methods give better accuracy than the explicit ones. For stiff problems, it is better to use Jacobian for solving the differential equations.

Matlab's NN Toolbox

Following program demonstrates the use of a Feed-forward network in Matlab

%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Matlab Example 1
% Learning XOR function using a MLP

% Create a FFN
net = newff([0 1; 0 1], [2 1], {'logsig', 'logsig'},'traingd');

%Training patterns
P = [0 0 1 1; 0 1 0 1];
T = [0 1 1 0];

% Training parameters
net.trainParam.show = 100;
net.trainParam.lr = 0.8;
net.trainParam.epochs = 3000;
net.trainParam.goal = 1e5;

%Train the network
[net,tr] = train(net, P, T);

% Test the network
Y = sim(net, P)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Following program demonstrates the use of RBFN in matlab to learn a nonlinear function

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MATLAB EXAMPLE 2
% Radial Basis Function Network (RBFN)
clear all;

% Data set for training

x = [0:4:200];
yd = sqrt(x);
spread = 20;

% Create a RBF network

net = newrb(x, yd, 1e-5,20);

%test data
P = [0:0.01:200];
T = sqrt(P);
T1 = sim(net, P);

%mean square error
mse(T1-T)
figure(2)
plot(P,T,'r', P, T1, 'b--');
legend('Desired output', 'Actual Output');
xlabel('Input');
ylabel('Output');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Matlab Editor : Rapid programming

Matlab has introduced something called "cells" to speed up writing/testing programs. Just go through this Flash animation to know more about cells.

http://www.mathworks.com/support/2007a/matlab/7.4/demos/RapidCodeIterationUsingCells_viewlet_swf.html


I had been avoiding Matlab for a long time. Its being a proprietory software. But I must admit, its simply great. They just need to reduce the price.

CAS: Roots of a polynomial

Maxima
---------------------------
(%i1) eq1: x^5-x+x*(1-x)^4*(5*a)+x^2*(1-x)^3*(10*b)+
x^3*(1-x)^2*10*(1-b)+x^4*(1-x)*5*(1-a)=0;
(%o1)

(%i2) solve(eq1,x);
(%o2)

(%i3) f1:rhs(%[1] );
(%o3)

(%i4) radcan(f1);
(%o4)

(%i5) ratexpand(%);
(%o5)

(%i6) d1:part(%,1);
(%o6)

(%i7) d2:part(%o5,[2,3,4]);
(%o7)

(%i8) ratsimp(%);
(%o8)

Roots in simplified format is given by

x = 0.5 +/- sqrt((b+1.5a-0.7)/(4b-2a-1.2))


Matlab
-----------------
syms a b x
f = x^5-x+x*(1-x)^4*(5*a)+x^2*(1-x)^3*(10*b)+
x^3*(1-x)^2*10*(1-b)+x^4*(1-x)*5*(1-a);
soln = solve(f,x);
f4 = soln(4,1);
f5 = soln(5,1);

f4 =
1/2/(5*a-10*b+3)*(5*a-10*b+3+(-75*a^2+
100*a*b-10*a+100*b^2-100*b+21)^(1/2))

f5 =
1/2/(5*a-10*b+3)*(5*a-10*b+3-(-75*a^2+
100*a*b-10*a+100*b^2-100*b+21)^(1/2))

It does not simplify roots any further. I tried all sorts of other commands like 'simplify', 'factor' etc. Is it possible to get the simplified form of roots as shown above?

Now scilab comes with its own symbolic toolbox and interface with 'maxima'. Check following link for details.
http://www.cert.fr/dcsd/idco/perso/Magni/s_sym/functions.html