Exercises

The beginning is simple

using Flux
using CUDA

A = randn(100,100)
b = randn(100)
A_g = A |> gpu
b_g = b |> gpu

To test the time, we measure the time for multiplication

julia> @time A*b;
0.069785 seconds (294.76 k allocations: 15.585 MiB, 14.75% gc time)

julia> @time A_g*b_g;
0.806913 seconds (419.70 k allocations: 22.046 MiB)

julia> @time A_g*b;
0.709140 seconds (720.01 k allocations: 34.860 MiB, 1.53% gc time)

We see that all three times are different. Can we infer anything from it? No! The problem is that during the first call to a function, some compilation usually takes place. We should always compare only the second time.

julia> @time A*b;
0.000083 seconds (1 allocation: 896 bytes)

julia> @time A_g*b_g;
0.000154 seconds (11 allocations: 272 bytes)

julia> @time A_g*b;
0.475280 seconds (10.20 k allocations: 957.125 KiB)

We conclude that while the computation on CPU and GPU takes approximately the same time, it takes much longer when using the mixed types.

To compare the results, the first idea would be to run

norm(A*b - A_g*b_g)

which would result in an error. We cannot use any operations on arrays stored both on CPU and GPU. The correct way is to move the GPU array to CPU and only then to compute the norm

julia> using LinearAlgebra

julia> norm(A*b - cpu(A_g*b_g))
1.2004562847861718e-5

The norm is surprisingly large. Checking the types

julia> (typeof(A), typeof(A_g))
(Matrix{Float64}, CUDA.CuMatrix{Float32})

we realize that one of the arrays is stored in Float64 while the second one in Float32. Due to the different number of saved digits, the multiplication results in this error.

We load the data, model and convert everything to GPU

using CUDA

m = Chain(
    Conv((2,2), 1=>16, relu),
    MaxPool((2,2)),
    Conv((2,2), 16=>8, relu),
    MaxPool((2,2)),
    flatten,
    Dense(288, size(y_train,1)),
    softmax,
)

file_name = joinpath("data", "mnist.bson")
train_or_load!(file_name, m)

m_g = m |> gpu
X_test_g = X_test |> gpu

Now we can measure the evaluation time. Remember that we need to compile all the functions by evaluating at least one sample before doing so.

m(X_test[:,:,:,1:1])
m_g(X_test_g[:,:,:,1:1])

julia> @time m(X_test);
1.190033 seconds (40.24 k allocations: 1.069 GiB, 21.73% gc time)

julia> @time m_g(X_test_g);
0.071805 seconds (789 allocations: 27.641 KiB)

Using GPU speeded the computation by more than ten times.

First, we load the data as many times before

m = Chain(
    Conv((2,2), 1=>16, relu),
    MaxPool((2,2)),
    Conv((2,2), 16=>8, relu),
    MaxPool((2,2)),
    flatten,
    Dense(288, size(y_train,1)),
    softmax,
)

file_name = joinpath("data", "mnist.bson")
train_or_load!(file_name, m)

When creating a table, we specify that its entries are Int. We save the predictions y_hat and labels y. Since we do not use the second argument to onecold, the entries of y_hat and y are between 1 and 10. Then we run a for loop over all misclassified samples and add to the error counts.

y_hat = onecold(m(X_test))
y = onecold(y_test)

errors = zeros(Int, 10, 10)
for i in findall(y_hat .!= y)
    errors[y[i], y_hat[i]] += 1
end

To create the dataframe, we use df = DataFrame(errors). It prints correctly integers and not strings. We change labels x1 to miss0, ... Similarly, we add the labels as the first column.

using DataFrames

df = DataFrame(errors, :auto)

rename!(df, [Symbol("miss$(i)") for i in 0:9])
insertcols!(df, 1, :label => string.(0:9))

To plot all these misclassified images, we find their indices and use the function imageplot. Since y are stored in the 1:10 format, we need to specify classes.

using ImageInspector

classes = 0:9

targets = onecold(y_test, classes)
predicts = onecold(m(X_test), classes)

imageplot(1 .- X_test, findall((targets .== 9) .& (predicts .== 7)); nrows=3)
savefig("miss.svg") # hide

To create the network and to load the data, we use

m = Chain(
    Conv((2,2), 1=>16, sigmoid),
    MaxPool((2,2)),
    Conv((2,2), 16=>8, sigmoid),
    MaxPool((2,2)),
    flatten,
    Dense(288, size(y_train,1)),
    softmax,
)

file_name = joinpath("data", "mnist_sigmoid.bson")
train_or_load!(file_name, m)

Before plotting, we perform a for loop over the digits. Then onecold(y_train, classes) .== i creates a BitArray with ones if the condition is satisfied, and zeros if the condition is not satisfied. Then findall(???) selects all ones, and ???[1:5] finds the first five indices. Since we need to plot the original image, and the images after the second and fourth layer (there is always a convolutional layer before the pooling layer), we save these values into z1, z2 and z3. Then we need to access to desired channels and plot then via the ImageInspector package.

using ImageInspector

classes = 0:9
plts = []
for i in classes
    jj = 1:5
    ii = findall(onecold(y_train, classes) .== i)[jj]

    z1 = X_train[:,:,:,ii]
    z2 = m[1:2](X_train[:,:,:,ii])
    z3 = m[1:4](X_train[:,:,:,ii])

    kwargs = (nrows = 1, size = (600, 140))
    plot(
        imageplot(1 .- z1[:, :, 1, :], jj; kwargs...),
        imageplot(1 .- z2[:, :, 1, :], jj; kwargs...),
        imageplot(1 .- z2[:, :, end, :], jj; kwargs...),
        imageplot(1 .- z3[:, :, 1, :], jj; kwargs...),
        imageplot(1 .- z3[:, :, end, :], jj; kwargs...);
        layout = (5,1),
        size=(700,800)
    )
    savefig("Layers_$(i).svg")
end

We plot and comment on three selected digits below.