Lecture 12: PLU Decomposition

For MATH3976 students: in the assignment, we will be using a bitstype. The following creates a new type of precisely 128 bits, that is a subtype of AbstractFloat:



In [39]:

    
bitstype 128 Float128 <: AbstractFloat

To create a Float128, we need to reinterpret another 128-bit type. The easiest case is UInt128:



In [43]:

    
u_int=rand(UInt128)
f=reinterpret(Float128,u_int);

typeof(f)









    Out[43]:





Float128

We can manipulate f by reinterpreting it back to a UInt128:



In [44]:

    
reinterpret(UInt128,f)









    Out[44]:





0x4488af451cad267a956f4096d9ef0c36

We see that it has exactly 128 bits:



In [45]:

    
bits(u_int)









    Out[45]:





"01000100100010001010111101000101000111001010110100100110011110101001010101101111010000001001011011011001111011110000110000110110"

We will need to access subsequences of the bits. In the following example, we decompose a 32-bit unsigned integer into two 16-bit unsigned integers.



In [46]:

    
x=rand(UInt32)

bits(x)









    Out[46]:





"01011101000011011101100100011001"

The syntax x % UInt16 drops the first 16 bits, and returns the last 16 bits as a UInt16:



In [47]:

    
x_16 = x % UInt16   # drops the first 15 bits, and keep the last 16 bits

bits(x_16)  # same as the last 16 bits of x









    Out[47]:





"1101100100011001"

To get at the first 16 bits, we will need to shift the bits right. This is equivalent to dividing by two and dropping the extra bits:



In [48]:

    
bits(UInt32(div(x,2)))









    Out[48]:





"00101110100001101110110010001100"

But it is more convenient to use x >> k, which shifts the bits right by k:



In [49]:

    
x_shift = x >> 1  # shifts the bits of x by 1, dropping the rightmost bit
bits(x_shift)









    Out[49]:





"00101110100001101110110010001100"



In [50]:

    
x_shift = x >> 2  # shifts the bits of x by 2, dropping the rightmost two bits
bits(x_shift)









    Out[50]:





"00010111010000110111011001000110"

We thus get the first and last 16 bits as follows:



In [51]:

    
x_shift = x >> 16 
x_first = x_shift % UInt16
x_last = x % UInt16

bits(x)









    Out[51]:





"01011101000011011101100100011001"



In [52]:

    
bits(x_first) * bits(x_last)









    Out[52]:





"01011101000011011101100100011001"



In [54]:

    
bits(x_first),bits(x_last)









    Out[54]:





("0101110100001101","1101100100011001")

Matrix norms

Just like vectors, matrices have norms that measure their "length". The simplest example is the Fr\"obenius norm, defined for an $n \times m$ real matrix $A$ as

$$\|A\|_F = \|{\rm vec}(A)\|_2 = \sqrt{\sum_{k=1}^n \sum_{j=1}^m A_{kj}^2}$$

This is using Julia's vec notation, which converts a matrix to a vector:



In [58]:

    
vec([1 2 3; 4 5 6; 7 8 9])









    Out[58]:





9-element Array{Int64,1}:
 1
 4
 7
 2
 5
 8
 3
 6
 9

While this is the simplest norm, it is not the most useful. In a lecture we will describe which norm is used in Julia.

The important thing for us is that if $\|A\| = 0$ then $A = 0$, so it can be used to test if two matrices are equal: if $\|A-B\| \approx 0$, then $A \approx B$:



In [59]:

    
A=rand(5,5)

Q,R=qr(A)

norm(Q*R-A)









    Out[59]:





5.507977739769666e-16