Search code examples
julialinear-algebramathematical-optimizationlinear-programmingflops

Mathematical flop count of column based back substitution function ( Julia )


I am new to Linear Algebra and learning about triangular systems implemented in Julia lang. I have a col_bs() function I will show here that I need to do a mathematical flop count of. It doesn't have to be super technical this is for learning purposes. I tried to break the function down into it's inner i loop and outer j loop. In between is a count of each FLOP , which I assume is useless since the constants are usually dropped anyway.

I also know the answer should be N^2 since its a reversed version of the forward substitution algorithm which is N^2 flops. I tried my best to derive this N^2 count but when I tried I ended up with a weird Nj count. I will try to provide all work I have done! Thank you to anyone who helps.

function col_bs(U, b)


n = length(b)
x = copy(b)

for j = n:-1:2
    if U[j,j] == 0
        error("Error: Matrix U is singular.")
    end
        x[j] = x[j]/U[j,j] 
        
        for i=1:j-1
        x[i] = x[i] - x[j] * U[i , j ]
        end
end

x[1] = x[1]/U[1,1]
 

return x
end

1: To start 2 flops for the addition and multiplication x[i] - x[j] * U[i , j ]

The $i$ loop does: $$ \sum_{i=1}^{j-1} 2$$

2: 1 flop for the division $$ x[j]  / = U[j,j] $$
3: Inside the for $j$ loop in total does: $$ 1 + \sum_{i=1}^{j-1} 2$$
4:The $j$ loop itself does:$$\sum_{j=2}^n ( 1 + \sum_{i=1}^{j-1} 2)) $$
5: Then one final flop for $$  x[1] = x[1]/U[1,1].$$

6: Finally we have 
$$\\ 1 + (\sum_{j=2}^n ( 1 + \sum_{i=1}^{j-1} 2))) .$$

Which we can now break down.

If we distribute and simplify
$$\\ 1 + (\sum_{j=2}^n + \sum_{j=2}^n \sum_{i=1}^{j-1} 2) .$$

We can look at only the significant variables and ignore constants,

$$\\
  \\ 1 + (n + n(j-1)) 
  \\ n + nj - n
  \\ nj
$$

Which then means that if we ignore constants the highest possibility of flops for this formula would be $n$ ( which may be a hint to whats wrong with my function since it should be $n^2$ just like the rest of our triangular systems I believe)

Function picture

Proof picture 1

Proof picture 2 and conclusion


Solution

  • Reduce your code to this form:

    for j = n:-1:2
       ...
       for i = 1:j-1
          ... do k FLOPs
       end
    end
    

    The inner loop takes k*(j-1) flops. The cost of the outer loop is thus

    \sum_{j=2}^n k (j-1)

    Since you know that j <= n, you know that this sum is less than (n-1)^2 which is enough for big O.

    In fact, however, you should also be able to figure out that

    \sum_{j=1}^n j = n (n+1) / 2