Solutions 7 - Examples in Sparse + Low-Rank Splitting

Assignment 1

We need functions from the notebook L13 Sparse + Low-Rank Splitting.


In [1]:
using LinearAlgebra
# Shrinkage
function Shr(x::Array{T},τ::T) where T
    sign.(x).*max.(abs.(x).-τ,zero(T))
end

# Singular value thresholding
function D(A::Array{T},τ::T) where T
    # This can be replaced by a faster approach
    V=svd(A)
    S=Shr(V.S,τ)
    k=count(!iszero,S)
    return (V.U[:,1:k]*Diagonal(S[1:k]))*V.Vt[1:k,:]
end


Out[1]:
D (generic function with 1 method)

In [2]:
function PCPAD(A::Array{T}) where T
    # Initialize
    δ=1.0e-7
    tol=δ*norm(A)
    m,n=size(A)
    S=zero(A)
    Y=zero(A)
    L=zero(A)
    T₁=zero(A)
    μ=(m*n)/(4*(norm(A[:],1)))
    μ₁=one(T)/μ
    λ=one(T)/sqrt(max(map(T,m),n))
    λμ₁=λ*μ₁
    ν=1e20
    maxiter=1000
    iterations=0
    # Iterate
    while (ν>tol) && iterations<maxiter
        # println(iterations," ",ν)
        iterations+=1
        L.=D(A-S+μ₁.*Y,μ₁)
        S.=Shr(A-L+μ₁.*Y,λμ₁)
        T₁.=(A-L-S)
        Y.+=(μ.*T₁)
        ν=norm(T₁)
    end
    L,S, iterations
end


Out[2]:
PCPAD (generic function with 1 method)

In [3]:
# For compilation
A0=rand(3,3)


Out[3]:
3×3 Array{Float64,2}:
 0.850396  0.232805   0.741312
 0.519471  0.0900631  0.393503
 0.490842  0.484163   0.517704

In [4]:
L,S,iter=PCPAD(A0)


Out[4]:
([0.7028471974285468 0.23280524268954814 0.7413119449568926; 0.37308537342236686 0.12357768690163247 0.3935032320233737; 0.49084212960961326 0.16258245252181958 0.5177044670576171], [0.14754859839422618 0.0 0.0; 0.14638525273159897 -0.03351459771798471 0.0; 0.0 0.3215807940816613 0.0], 99)

In [5]:
rank(L), L+S-A0


Out[5]:
(1, [0.0 -7.996942155696907e-8 4.9466783336171716e-8; 0.0 -4.163336342344337e-17 -1.0293989405329995e-7; 6.280581277273711e-8 -2.220446049250313e-16 -2.2554651368800194e-8])

In [6]:
L


Out[6]:
3×3 Array{Float64,2}:
 0.702847  0.232805  0.741312
 0.373085  0.123578  0.393503
 0.490842  0.162582  0.517704

In [7]:
# Packages for video manipulation
# import Pkg; Pkg.add("VideoIO")
# import Pkg; Pkg.add("Makie")

In [7]:
using Images
using Makie, VideoIO

In [8]:
varinfo(VideoIO)


Out[8]:
name size summary
VideoIO 1.015 MiB Module
appendencode! 0 bytes typeof(appendencode!)
encode! 0 bytes typeof(encode!)
encodevideo 0 bytes typeof(encodevideo)
finishencode! 0 bytes typeof(finishencode!)
gettime 0 bytes typeof(gettime)
mux 0 bytes typeof(mux)
opencamera 0 bytes typeof(opencamera)
openvideo 0 bytes typeof(openvideo)
play 0 bytes typeof(play)
playvideo 0 bytes typeof(playvideo)
prepareencoder 0 bytes typeof(prepareencoder)
pump 0 bytes typeof(pump)
read 0 bytes typeof(read)
read! 0 bytes typeof(read!)
viewcam 0 bytes typeof(viewcam)

In [9]:
video=VideoIO.open("files/visor.avi")


Out[9]:
AVInput(files/visor.avi, ...), with
  1 video stream(s)

In [10]:
playvideo(video)

In [11]:
VideoIO.get_duration("files/visor.avi")


Out[11]:
229.2

In [12]:
# Read the frames in Gray
f = VideoIO.openvideo("files/visor.avi",target_format=VideoIO.AV_PIX_FMT_GRAY8)
frames=Array{Array{Gray{Normed{UInt8,8}},2},1}(undef,0)
while !eof(f)
    push!(frames,read(f))
end
close(f)

In [13]:
# Number of frames
length(frames)


Out[13]:
2292

In [14]:
# See particular frame
frames[1401]


Out[14]:

In [15]:
# Make shorter video clip
clip=frames[1401:1450]


Out[15]:
(a vector displayed as a row to save space)

In [16]:
size(clip[1])


Out[16]:
(240, 320)

In [17]:
# Turn frames into tall matrix
mi,ni=size(clip[1])
m=mi*ni
n=length(clip)
A=Array{Float64}(undef,m,n)
for i=1:n
    A[:,i]=vec(float(clip[i]))
end

In [18]:
size(A)


Out[18]:
(76800, 50)

In [19]:
norm(A)


Out[19]:
993.2119585880553

In [21]:
# For orientation
@time Q=qr(A);
@time D(A,0.5);


  0.112519 seconds (5 allocations: 29.325 MiB)
  0.431027 seconds (25 allocations: 146.607 MiB, 4.50% gc time)

In [23]:
# Compute the splitting - 4 minutes
@time L,S,iters=PCPAD(A)


535.896983 seconds (42.02 k allocations: 361.787 GiB, 7.74% gc time)
Out[23]:
([0.18823530800601174 0.18823530357649265 … 0.17647056764610916 0.17647054229722997; 0.18823530800602248 0.18823530357650342 … 0.17647056764611407 0.1764705422972349; … ; 0.7176470909901755 0.7176470920849694 … 0.7058823769950273 0.705882400571486; 0.7019607808997179 0.7019607832012424 … 0.7137255012368067 0.7137254736428041], [0.0 0.0 … -0.0 -0.0; 0.0 0.0 … -0.0 -0.0; … ; 0.0 0.0 … -0.0 -0.0; -0.0 -0.0 … -0.0 -0.0], 1000)

In [24]:
# Reconstruct the low-rank video component
rank(L), norm(A-L-S), norm(A)


Out[24]:
(16, 0.0004293692397533305, 993.2119585880553)

In [25]:
svdvals(L)[1:20]


Out[25]:
20-element Array{Float64,1}:
 1000.685641671293
   14.872151740792042
   11.379472221987633
    9.729815137981031
    7.706656406174793
    4.156856638710673
    3.180028276068219
    2.45653107143592
    1.2897504608196235
    0.7788793555613616
    0.6411961879091401
    0.5972429272521264
    0.01857625930341237
    0.008823112225258835
    0.000563932220783042
    0.0002809520311917548
    1.4467087136026894e-13
    2.7156280425745636e-14
    2.7056621945686743e-14
    2.7016194859395314e-14

In [26]:
# How to restore the video?
# First frame of the low-rank part
v1=L[:,1]
v2=reshape(v1,mi,ni)
map(Gray{Normed{UInt8,8}},v2)


Out[26]:

In [27]:
# First frame of the sparse part 
s1=S[:,1]
s2=reshape(s1,mi,ni)
map(Gray{Normed{UInt8,8}},clamp01.(s2.+0.5))


Out[27]:

In [28]:
LowRank=similar(clip)
Sparse=similar(clip)
for i=1:n
    LowRank[i]=reshape(L[:,i],mi,ni)
    Sparse[i]=clamp01.(reshape(S[:,i],mi,ni).+0.5)
end

In [29]:
# import Pkg; Pkg.add("JLD")

In [22]:
# Let us save the results
using JLD

In [31]:
@save "files/visor_results.jld" LowRank Sparse

In [23]:
@load "files/visor_results.jld"


Out[23]:
7-element Array{Symbol,1}:
 :LowRank
 :Sparse
 Symbol("_creator\\ENDIAN_BOM")
 Symbol("_creator\\JULIA_MAJOR")
 Symbol("_creator\\JULIA_MINOR")
 Symbol("_creator\\JULIA_PATCH")
 Symbol("_creator\\WORD_SIZE")

In [24]:
LowRank[1]


Out[24]:

In [25]:
Sparse[1]


Out[25]:

In [26]:
# Play the LowRank part a bit slower (framerate=10)
props = [:priv_data => ("crf"=>"22","preset"=>"medium")]
encodevideo("files/LowRank.mp4",LowRank,framerate=10,AVCodecContextProperties=props)


┌ Info: Video file saved: C:\Users\Ivan_Slapnicar\Documents\GitHub\GIAN-Applied-NLA-Course\src\Module C - Applications/files/LowRank.mp4
└ @ VideoIO C:\Users\Ivan_Slapnicar\.julia\packages\VideoIO\TAILL\src\encoding.jl:221
┌ Info: frame=   50 fps=0.0 q=-1.0 Lsize=      54kB time=00:00:04.70 bitrate=  94.1kbits/s speed=2.36e+03x    
└ @ VideoIO C:\Users\Ivan_Slapnicar\.julia\packages\VideoIO\TAILL\src\encoding.jl:222
┌ Info: video:53kB audio:0kB subtitle:0kB other streams:0kB global headers:0kB muxing overhead: 1.879838%
└ @ VideoIO C:\Users\Ivan_Slapnicar\.julia\packages\VideoIO\TAILL\src\encoding.jl:223
Out[26]:
"files/LowRank.mp4"

In [28]:
videoLowRank=VideoIO.open("files/LowRank.mp4")
playvideo(videoLowRank)

In [29]:
# Play the Sparse part
encodevideo("files/Sparse.mp4",Sparse,framerate=10,AVCodecContextProperties=props)


┌ Info: Video file saved: C:\Users\Ivan_Slapnicar\Documents\GitHub\GIAN-Applied-NLA-Course\src\Module C - Applications/files/Sparse.mp4
└ @ VideoIO C:\Users\Ivan_Slapnicar\.julia\packages\VideoIO\TAILL\src\encoding.jl:221
┌ Info: frame=   50 fps=0.0 q=-1.0 Lsize=      97kB time=00:00:04.70 bitrate= 169.8kbits/s speed=N/A    
└ @ VideoIO C:\Users\Ivan_Slapnicar\.julia\packages\VideoIO\TAILL\src\encoding.jl:222
┌ Info: video:96kB audio:0kB subtitle:0kB other streams:0kB global headers:0kB muxing overhead: 1.032880%
└ @ VideoIO C:\Users\Ivan_Slapnicar\.julia\packages\VideoIO\TAILL\src\encoding.jl:223
Out[29]:
"files/Sparse.mp4"

In [31]:
videoSparse=VideoIO.open("files/Sparse.mp4")
playvideo(videoSparse)

In [ ]: