[TECH] Maintaining a Load Factor (α) in Dynamic Data structures with constant (Θ(1)) amortized cost per operation.

Accounting Method is an elegant technique for amortized analysis of the data structures. Most of us should have encountered Dynamic data structures (which need to expand and contract to maintain some invariant) at some point of time. A stack implementation is one example of a dynamic data structure. As we perform pop operation on the stack eventually it become full (over-flow), at this stage we might choose to increase (may be by doubling or some other heuristic) the size of the stack. We also want to handle the under-flows i.e. if the number of elements in the stack is too small we might want to relocate the elements into a new stack ( with smaller size ). In general the dynamic data structures have a parameter called the load factor and our aim is to maintain that. One simple way to maintain this is to handle under-flows and over-flows on demand. Unfortunately this will not give us a amortized cost per operation on the stack. However we can still get amortized cost by avoiding this on demand expansion and compaction.

Let denote the number of elements present on the stack after operations (including push and pop). Let denote the allocated size of the stack after operations. We Define a potential energy function as follows.

Clearly . The amortized cost of the push operation which does not require an expansion of the stack is as follows.

The amortized cost of the push operation which requires doubling of the array is as follows.

The amortized cost of the pop

From the above we can see that that amortized cost of all the operations is when we use the doubling heuristic for the stack.

(b)

Let the load factor and let . We use the following heuristic to handle expansion and compaction of the array.

"When (i.e. the array is full) increase the size of the array by times. When (i.e. the array about to under flow) reduce the size of the array by half ()".

We define the following potential function and finally show that the amortized complexity of each of the operations push, pop is .

We now prove that the amortized cost of each operation is bounded by a constant. From the previous problem (since the potential function is same as the previous problem when ) we can observe that push operation when (expansion) has an amortized cost of . We can also observe that whenever expansion does not occur the amortized cost of the push operation is . So we will concentrate on the other possible situations. We use the notation to indicate the load factor of the array after operations on the array.

Amortized cost of push when

Amortized cost of push when ,

Amortized cost of pop when

Amortized cost of pop when

Amortized cost of pop when and

Amortized cost of pop when with compaction

If the operation is a pop and we find that the load factor then we apply compaction. Compaction reduces the size of the allocated array by . So after compaction . The amortized cost of the pop and compaction is as follows.

We have proved that the amortized cost of both push and pop with all possible values of the load factor to be .

[TECH] Exponential function (e-μt) is the only unique solution satifying R(x+t) = R(x)R(t)(semi-group property).

Exponential function keeps popping up every where stochastic analysis. However it has a crucial property which is often exploited during the analysis. The property I'm referring to is R(X+Y) = R(Y)R(Y) (e.g. X and Y) could be random variables). Now I prove a strong statement " R(X+Y) = R(X)R(Y) If and Only If R(x) = e^-μt where μ is some constant. The reverse direction (<==) is trivial because we can easily verify the fact by substituting R(x) with e^-μt. In the following we prove that other direction (==>).

One may wonder how about function a^x which also satisfies the semi-group property by inspection. However the μ takes care of that observe that a^x = e^ln(a)x.

[TECH] Simple comparision of real numbers

We know that the fixed precision representation of real numbers is done by mapping them into integer space I from real space R. However the underlying integer representation of the 32-bit floats and 64-bit doubles is a sign magnitude representation (that is the reason why we see -0.00 during output of some numerical algorithms). People seem to be adding a constant to convert the underlying sign magnitude form to 2's compliment. However following is a more logical conversion (since we know that 2's compliment is 1's compliment +1 , the motivation of 2's compliment is to avoid the -0 as we know). The following comparison routine compares two floating point numbers (add the tolerance in the terms of the number the real numbers to make it approximate comparison).

/*Less than or greater than*/
char CompareFloats(float a, float b){
    int aint = *(int *)&a;
    int bint = *(int *)&b;
    aint = (aint & (1UL<<31))? ~(aint^(1UL<<31))+1:aint;
    bint = (bint & (1UL<<31))? ~(bint^(1UL<<31))+1:bint;
    aint -= bint;
    printf("%f %s %f \n",a,(aint <0)?"<":">=",b);
    printf("There are %d real numbers \n",(aint<0)?-aint:aint);
    printf("between %f and %f (in 32 bit float representation)\n",a,b);
}

[TECH] Launchpad and Bazaar

I have now started to move all my development efforts to launchpad and bazaar. I really found launchpad a great place to maintain and release my work. I tried several alternatives before I moved to launchpad, I tried to use code.google but unfortunately its not as rich as launchpad and bazaar. I have released version 0.9 of libEditScript project whose aim to build a high performance and space efficient sequence alignment library. I have charted out the features for the next release, one thing I'm sure people would like to use is a JNI (Java Native Interface) so that people can use this space efficient code in their java code. I have seen that this problem of the need of edit script arises in several occasions, I have seen people using the space inefficient version O(n^2) of edit script computation. This is were we gain significantly in terms of the space saving. Also I have one more idea to make it more parallel in the sense I want to use some concurrent techniques to make this parallel, especially since I have circular queue I can easily use several threads to make this concurrent which I will add as a next step in the next release..

Checkout the libEditScript at https://launchpad.net/libeditscript. Next in my list would be to make a release of the weightedmatching library.

[TECH] A non-recursive algorithm to compute the Edit Script in O(min(n1,n2)+log(min(n1,n2)) space

Recently I was working on an area efficient and synthesizable design to compute the edit script (minimum cost sequence of INSERT, DELETE and CHANGE) operations to transform string S₁ to S₂. After a little bit of thought I have the following idea to build a non-recursive version of the Hirschberg's algorithm, since the algorithm is non-recursive we can build an efficient digital circuit with this idea. We use a simple circular queue and apply DFS (Depth First Search) and we can prove that the capacity of this queue at any stage of the algorithm is θ(log(min(n₁,n₂)). The proof is simple we choose the geometrically decreasing string which has smaller length from the given strings(min(n₁,n₂)). Since we do a depth first search and the depth of subproblem tree is θ(log(min(n₁,n₂)), so we will have atmost θ(log(min(n₁,n₂)) subproblems in the circular queue at any stage of the algorithm.

The core non-recursive algorithm starts off with the initial problem and finds the minimum alignment split ('q') between (S₁[1.....n₁/2] and S₂ creates two sub problems and appends in the circular queue (currently it uses BFS). Download the linear space implementation from here Following is the outline of the core non recursive implementation.

     91     CQueue bfs_list;
     92     ESSubProblem es_sub;
     93 
     94     InitializeCQ(&bfs_list,sizeof(ESSubProblem));
     95     /*Put the toplevel subproblem into the CQueue*/
     96     es_sub.sub_s1 = s1; es_sub.sub_n1 = n1;
     97     es_sub.sub_s2 = s2; es_sub.sub_n2 = n2;
     98 
     99     AppendCQ(&bfs_list,&es_sub);
    100     while(DequeueCQ(&bfs_list,&es_sub)){
    101         sub_s1 = es_sub.sub_s1; sub_s2 = es_sub.sub_s2;
    102         sub_n1 = es_sub.sub_n1; sub_n2 = es_sub.sub_n2;
    103 
    104         if(sub_n1 <=1){
    105             /*If sub_n1 is 1 or 0 we cannot split it any more*/
    106             EditDistanceLS(sub_s1,sub_n1,sub_s2,sub_n2);
    107             i = sub_n1; j = sub_n2;
    108             /*Figure out the edit operations*/
    109             while(i || j){
    110                 op = GetOp(i,j,(i?sub_s1[i-1]:'\0'),
    111                     (j?sub_s2[j-1]:'\0'));
    112                 switch(op){
    113                     case 'I':
    114                             EditScriptS2[(&sub_s2[j-1]) - s2] = 'I';
    115                             j--;
    116                             break;
    117                     case 'D':
    118                             EditScriptS1[(&sub_s1[i-1]) - s1] = 'D';
    119                             i--;
    120                             break;
    121                     case 'C':
    122                             EditScriptS1[(&sub_s1[i-1]) - s1] =
    123                                 (sub_s1[i-1] == sub_s2[j-1])?':':'C';
    124                                 i--; j--;
    125                 }
    126             }
    127         }else {
    128             /*Add the two subproblems*/
    129             q = FindMinSplit(sub_s1,sub_n1,sub_s2,sub_n2);
    130             /*SUB PROBLEM 1*/
    131             es_sub.sub_n1 = CELING(sub_n1,2);
    132             es_sub.sub_n2 = q;
    133             AppendCQ(&bfs_list,&es_sub);
    134 
    135             /*SUB PROBLEM 2*/
    136             es_sub.sub_s1 = &(sub_s1[CELING(sub_n1,2)]);
    137             es_sub.sub_s2 = &(sub_s2[q]);
    138             es_sub.sub_n1 = sub_n1/2;
    139             es_sub.sub_n2 = sub_n2-q;
    140             AppendCQ(&bfs_list,&es_sub);
    141         }
    142     }

[TECH] A Graph with θ(log(n)) approximation ratio for greedy vertex cover

We know that the greedy algorithm (which picks the vertex with next highest degree) to solve the vertex cover problem is not optimal. However we have 2-Approximation algorithm for the vertex cover problem (based on a integer linear programming formulation). Actually it was not very straight forward to for me prove that the greedy algorithm has an approximation ratio greater than 2. However after some thinking the following construction indeed proves that the greedy algorithms approximation ratio is > 2. In fact we can easily see (from the figure) that the greedy algorithm will first pick the node with highest degree i.e n (the magenta). Then it deletes all the edges adjacent on the magenta node and continues and clearly the greedy algorithm will have all the nodes in the blue bounding box. However its clear that if we choose the nodes in the oval region we get a minimum vertex cover of size n. The size of greedy vertex cover is n/2 + n/3 + n/4 .... n/n ≈ θ(nlog(n)) (for large n) and hence the approximation ratio for this graph is θ(log(n)).

[TECH] On The Uniqueness Of a Perfect Match In Undirected Acyclic graphs.

I was studying about several results related to Tiling Theory. Tiling Theory adds formalism to the question " Given a simple rectangular region (R) can we find a valid tiling by using tiles from a set (T) ". If our tile set (T) is a singleton set with a unit square then we can tile any region. However what if the the tile set contains only a 2x1 rectangle (domino)? what if the tile set contains a set of polymino's (tetris game)?. Although the question may seem very simple there is a rich combinatorial work based on this question. In fact the study about the tiling problem originated from Wang's Conjecture about periodic tiling. Actually the algorithmic research about tiling theory was dormant for quite some time, even the proof for 2x1 dominoes is fairly recent see this . I was reading about the proof and was really impressed by its elegance. Often people consider simple things useless or don't have enough respect for simple facts, however I always find that good results are based on several very simple facts. In fact I would like to highlight simple ideas in every algorithm I study about.

In this proof to prove that the problem is NP-complete for 2x1 domino's a simple fact that " There is exactly only one perfect matching (if at all one exists) in a Tree " is used. This can be proved intuitively as follows.

Let M be the perfect match in the Tree T, if M is not a unique perfect match then there should be some other perfect match lets call it M^'. If we observe the leaf nodes in T which have only one edge adjacent on them so both M and M^' should agree on the matched edges which have leaf nodes of T as the vertices . We can now delete all the leaf nodes and their adjacent edges and continue the argument until we don't have any leaf nodes. Thus from this argument its clear that both M and M^' are the same and hence there is only a unique perfect match if at all if one exists in a tree T.

[TECH] Converting Market Matrix Sparse Representation to Row Compressed Sparse Representation

I wanted to use some good benchmarks to illustrate the performance gains of the Fast Dual Update Algorithm for pre-ordering sparse matrices. University of Florida Sparse Library has a rich collection of sparse matrices from various domains. Unfortunately all the matrices are in Market Matrix format but not in Row Compressed format which I currently use.

The following is a handy script which converts the Market Matrix format to Row Compressed Format mm2rc.pl. The script uses a simple logic based on external sorting.

[TECH] Computing Eigen Values and Eigen Vectors using LAPACK

I have used the LAPACK to compute the eigen values, till recently I needed the eigen vectors both left and right. I thought this could be really handy if you need to compute eigen values and eigen vectors. Download both library and this file from here

/*Computing both the eigen values and eigen vectors using LAPACK.
 *vamsi.krishnak (at) gmail (dot) com
 **/
#include "f2c.h"
#include "clapack.h"
#include<stdio.h>
#include<math.h>
#include<malloc.h>
#include<assert.h>

int main(){
    integer dim=0;
    unsigned int i;
    doublereal *A = NULL;
    real *eig_values = NULL;
    real *work;
    char JOBVL='V'; /*Left eigen vectors*/
    char JOBVR='V'; /*No right eigen vectors*/
    integer iwork[1];
    integer lwork,liwork=1,info,lda;
    doublereal *WR=NULL;
    doublereal *WI=NULL;
    doublereal *VL=NULL;
    doublereal *VR=NULL;
    integer LDVR=1;
    integer LDVL=1;
    integer LDA;
    doublereal *WORK=NULL;
    integer LWORK;


    printf("Please Enter the dimension of the array:\n");
    scanf("%d",&dim); LDA=dim; lda=dim; LDVL=dim;LDVR=dim;
    printf("Please enter the Matrix A\n");
    eig_values = (real *) malloc(sizeof(real)*dim);
    lwork = 10*dim; 
    work = (real *) malloc(sizeof(real)*lwork);
    A = (doublereal *) malloc(sizeof(doublereal)*(dim*dim));
    printf("Please Enter the elements of the matrix\n");
    i = 0;
    
    /*Double real and imaginary parts*/
    WR = (doublereal *)malloc(sizeof(doublereal)*dim);
    WI = (doublereal *)malloc(sizeof(doublereal)*dim);
    VL = (doublereal *)malloc(sizeof(doublereal)*(dim*dim));
    VR = (doublereal *)malloc(sizeof(doublereal)*(dim*dim));
    LWORK = 6*dim;
    WORK = (doublereal *)malloc(sizeof(doublereal)*LWORK);  
    assert(WI && WR && VL && VR);
    /**************/
    do{
        scanf("%lf",&A[i]);
        i++;
    }while(i<(dim*dim));
    /*Compute the Eigen values*/
    //ssyev_("N","U",&dim,A,&dim,eig_values,work,&lwork/*,iwork,&liworki*/,&info);
    dgeev_(&JOBVL,&JOBVR,&dim,A,&LDA,WR,WI,VL,&LDVL,VR,&LDVR,WORK,&LWORK,&info);
    if(!info){
        for(i=0;i<dim;i++){
            printf("%lf+i%lf\n",WR[i],WI[i]);
        }
        printf("===Right Eigen Vectors===\n");
        int k; int conj;
        for(i=0;i<dim;i++){
            /*Print i^th right eigen vector*/
            if(fabs(WI[i] > 0)){
                /*eigen vector i and eigen vector i+1*/
                for(conj=0;conj<2;conj++){
                    printf("[ ");
                    for(k=0;k<dim;k++){
                        if(!conj){
                            printf("(%lf)+i(%lf) ",VL[i*dim+k],VL[(i+1)*dim+k]);

                        }else{
                            printf("(%lf)-i(%lf) ",VL[i*dim+k],VL[(i+1)*dim+k]);
                        }
                    }
                    printf("  ]\n");
                }
                i++;
            }else{
                /*real eigen vector*/
                printf("[  ");
                for(k=0;k<dim;k++){
                    printf("%lf ",VL[i*dim+k]);
                }
                printf("] \n");
            }
            printf("\n");
        }
    }
}

[TECH] Cauchy's Inequality [continued...]

Just adding the solution to one more problem continuing from last post.

[TECH] Cauchy's inequality a great tool for approximation algorithms

Its been really long since I made a post, all these day's I'm doing a lot of theory and trying to come up with approximation algorithms to the Border Length Minimization Problem. After a long time I did prove the problem is NP-hard by reducing the T.S.P problem to the B.L.M.P problem. I'm now trying to find some approximation algorithms for the problem. I thought if I could strengthen my skills on inequalities it would be really great and I found a great problem book for inequalities by Michael Steele which has a good collections of problems on inequalities. In this post I provide my solutions to some of the problems which seemed interesting.