EMTH211 Assignment 2

Joshua SMIT: 91004964

Jacob OGILVIE: 34434903

Question 1a: Computing Power

The data in Table 1 gives the number of floating point operations per second (FLOPS) that
the world’s fastest supercomputer was capable of each year between 1993 and 2024. The
values span several orders of magnitude, so a logarithmic scale is used on the 𝑦-axis.

Figure 1: FLOPS of the world’s fastest supercomputer (1993–2024).

Question 1b:

We aim to fit a curve of the form:

𝑦 = 𝑎𝑒 !"
Where x is the year and y is the corresponding FLOP rate. Taking the natural logarithm
gives:
𝑙𝑛 𝑦 = 𝑙𝑛 𝑎 + 𝑏 𝑥
This is a linear relationship between 𝑙𝑛 𝑦 and x. Letting 𝑧 = 𝑙𝑛 𝑦, the model becomes linear
and can therefore be expressed in the matrix form:

𝑨𝒑 ≈ 𝒛

where,

𝟏 𝔁𝟏 𝒍𝒏 𝒚𝟏
𝟏 𝔁𝟐 𝒍𝒏 𝒂 𝒍𝒏 𝒚𝟐
𝑨 = 0 4 𝒑 = 5 : 𝒛 = 0 4
⋮ ⋮ 𝒃 ⋮
𝟏 𝔁𝒏 𝒍𝒏 𝒚𝒏
Here, the matrix 𝐴 has a column of ones and a column containing the year values. The vector
𝑝 contains the unknown values 𝑙𝑛 𝑎 and 𝑏. The vector 𝑧 contains the natural logarithms of the
measured flop rates.

We can use a QR decomposition to solve this least square problem. Where 𝑄 is a orthogonal
matrix and 𝑅 is upper triangular. Substituting this into the linear system gives,
𝑄𝑅𝑝 ≈ 𝑧
Because 𝑅 is upper triangular, this system can be solved for 𝑝 using back substitution. This
gives us the two unknows parameters contained in 𝑝. Once we have the numerical values of
𝑙𝑛 𝑎 and 𝑏, we can reconstruct the original exponential model.
𝑦 = 𝑒 ('()*!") = 𝑎𝑒 (!")
Finally, the fitted model is plotted with the measured data on a semi log y-axis, where it
appears as a straight line.

Figure 2: Exponential least squares fit plotted alongside the measured FLOPS data using a semilogarithmic y-axis.

Question 1c:

To improve the accuracy of the model, the dataset was divided into two segments at 2011,
this was chosen by visual inspection as a point where the growth trend appears to change.
Two different least squares fits were then applied to the early (𝑥 < 2011) and late ( 𝑥 > 2011)
segments. For each segment like above the data was linearized using,
𝑙𝑛 𝑦 = 𝑙𝑛 𝑎 + 𝑏 𝑥,
and least squares fit was computed using QR decomposition. This produced two separate
exponential models, each representing the growth behavior in their respective period.
The fitted curves and measured data were plotted together on a semi log y-axis.
Splitting the data in this way is justified because there is clear evidence of a change in the
growth rate around 2011. This period corresponds to the rise in GPU accelerated computing.
Developments in chip architecture also shifted towards multi core computed, which
drastically impacted the computers performance. As a result, modelling the two time periods
separately gives a better representation of the data.

Figure 2: Two-piece exponential least squares fit (split at 2011) illustrating different growth rates before and after the split.

Question 1d:

To evaluate the accuracy of our least squares models, the error between the actual and
predicted FLOP values was computed using the formula:
Error = |ln(𝑦),-.)' ) − ln(𝑦/0123,-12 )|,
where 𝑦),-.)' is the measured FLOPs for a given year and 𝑦/0123,-12 is the value predicted
by the model. For a single exponential fit the maximum error within the dataset occurs in
2013 with a value of 1.19, measured in logarithmic units. For the two-piece fit, this reduces
the maximum error to 0.62 occurring in 2005, which reflects a more accurate representation
of the data. Given a FLOPS value of 3.600x104 , using the single exponential fit, the
predicted FLOPS for 1969 is 1.84x105 , with an error of 5.28. Using the ‘Early’ segment of
the two line fit (𝑥 < 2011) the predicted FLOPS for 1969 is 1.53x106 , with an error of 7.77.

Both predictions are significantly less accurate than the errors observed within the fitted data
range. This highlights the limitations of extrapolation compared to interpolation; while the
models perform well for years within the dataset, extending them outside the observed range
leads to large deviations, on the order of several magnitudes because the underlying growth
trends are unlikely to have remained constant over such a long period.

Question 2a: Singular Value Decomposition (SVD)

For this question I adapted the given skeleton code to create the following 4 Figures. These
Figures use SVD to decompose A as:

𝐴 = 𝑈Σ𝑉 7

Where A is a matrix of brightness values for either the Red, Green or Blue colour channels for
each pixel of the given photo. The final image is a combination of these shaded values
combining to the colour that we see. These values are converted to a float, or decimal and this
would be visualized in Table 1 below.

RED GREEN BLUE Output

1 0 0 Pure Red
0 1 0 Pure Green
0 0 1 Pure Blue
1 1 0 Yellow
1 0 1 Purple
0 0 0 Black
1 1 1 White

Table 1: Colour Combinations

A is then decomposed down to U, å, and VT. U is the basis for the row patterns. This explains
how the variation of the row’s changes, and VT is the same for the columns. å then, tells us
how significant each of these patterns are and how much they contribute to the image. The
following listing shows how my code decomposes the original matrix.

def svd_compress_image(image, k):

""" Compresses the image using SVD """

comp = np.zeros_like(image)
for i in range(3):
U, S, Vt = [Link](image[..., i}, full_matrices=False)
comp[..., i] = (U[:, :k] @ [Link](S[:k]) @ Vt[:k, :])
return [Link](comp, 0, 1)

Listing 1: SVD Function

I used the two following images for decomposing with a k value (rank) of 100,

Figure 4: IMG_1682 Figure 5: IMG_1626

Using a Rank (100) approximation,

Figure 6: Rank (100) Approximation of given Images

As seen by the figures above, using the Rank (100) approximation we get almost an exact
recreation. This is because the code keeps the 100 largest singular values and their
corresponding patterns. This means that when the image is recreated by these metrics, the
image is very similar to the original, with slightly less fine detail.
Question 2b:
I again ran the same code but with different values for k, these being [1, 10, 50], and the
results are below.

Figure 7: SVD for k = 1

Figure 8: SVD for k = 10

 Figure 9: SVD for k = 50

After reviewing these images, I would say 10 singular values gives a vaguely recognizable
approximation of the image, whereas at 50, the image is a good representation of the original
image, again missing the finer detail, but with enough of a recognizable pattern to discern the
original image. As seen in the images, they seem to retain a similar definition by each different
rank. I believe this is due to the similarities of colour, and intricacy of the images. Both have
similar detail and sharp changes in patterns.

Question 2c:
For this part of the question, we are asked to find the how many real numbers are needed to
specify the first k singular values and vectors to construct the approximation, and how many
real numbers are needed to represent the array of the original image. From our code, we load
the image using the lines,

image = [Link](filename).astype(np.float64)
if [Link]() > 1.0:
image /= 255.0

Listing 2: Image loading

to and dividing by 255.0 ensures that the value within the array is a decimal between 1 – 0. The
shape of the image is given by,

[Link] = (height, width, 3)

Listing 3: Shape definition

which shows the array is defined by the height, width and the value 3, representing Red, Green
and Blue colour channels. This means for the original image the total number of real numbers
needed is given by,

𝑁8039:()' = 3 × ℎ𝑒𝑖𝑔ℎ𝑡 × 𝑤𝑖𝑑𝑡ℎ.

Equation 1: Real Numbers in Original Image

Running the following code shown in Listing 4, prints the number of rows and columns as well
as the total number.

image = load_image("IMG_1626.JPG")
rows, cols, channels = [Link]

print(f"Image dimensions: {rows} rows × {cols} columns")

print(f"Number of colour channels: {channels}")
print(f"Total number of real numbers: {rows * cols * channels}")

Listing 5: Rows and Columns

The printed output for the code in Listing 5, for each image is shown below in Table 2.

Rows Columns Channels NOrigonal

IMG_1626 3,024 4,032 3 36,578,304
IMG_1682 1,203 679 3 2,450,511
Table 2: Code Snippet Outputs

This however changes for the compressed images since we are finding the SVD. This means
our approximation is given by,

𝐴; = 𝑈; 𝑆; 𝑉;7 ,

Where k is the rank of the approximation. Due to this we need to store three things,

1. 𝑈; : an m x k matrix – mk numbers
2. 𝑉;7 : an n x k matrix – nk numbers
3. 𝑆; : a vector of k singular values – k numbers

Due to this the formula to find the amount of real numbers stored for our SVD approximations
of the images can be given by,

𝑁<=> = 3𝑘(𝑚 + 𝑛 + 1).

Equation 2: Real Numbers in SVD Approximation
From here we can also find the efficiency of each rank approximation. This is found by dividing
the amount of real numbers from each SVD by NOrigonal. The results of this equation for each
image and the corresponding k values are displayed in the tables below.

NSVD (k = 1) Efficiency (%)

IMG_1626 21,171 0.058
IMG_1682 5,649 0.231
NSVD (k = 10) Efficiency (%)
IMG_1626 211,710 0.578
IMG_1682 56,490 2.305
NSVD (k = 50) Efficiency (%)
IMG_1626 1,058,550 2.894
IMG_1682 282,450 11.52
NSVD (k = 100) Efficiency (%)
IMG_1626 2,117,100 5.788
IMG_1682 564,900 23.05
Table 3: Comparison of Real Numbers need per Rank

From these tables we can see that as k increases so does the efficiency, this means however,
that as we increase the definition, we increase the amount of numbers we need to store
compared to the original image. The trade-off is clear: the rank-100 approximation achieves
nearly identical visual quality to the original image while storing only a small fraction of the
total data. This occurs because most of the important image information is captured by the
largest singular values, which represent the dominant patterns and structures in the image,
while the smaller singular values contribute mainly to fine detail and noise.

Question 2d:
When compressing a matrix, each colour channel is treated as a matrix and compressed with
SVD. Relating this back to Eq. 2 in 2c, by changing the rank, k, for each colour channel we
can decrease the overall storage required for the image while keeping definition almost
consistent. This is further impacted by the idea that the human eye is more sensitive to green
than it is blue. Due to this the rank of the blue colour channel could be a lot less than the green
while retaining most, if not all, of the fine detail provided by the image. Below in Figure 10
and Figure 11 we see the comparison of two different images and how they differ with various
ranks for different colour channels. These reduced images are compared to the original but for
the Rank(50) approximation, IMG_1626 uses » 3% and IMG_1682 uses » 12% of the numbers
stored respectively by the original.
 Figure 10: Comparison of Compressed IMAGE_1626 with SVD

Figure 11: Comparison of Compressed IMAGE_1682 with SVD

As we can see by these figures, even though they are significantly reduced in storage size, they
produce vary similar images. MSE is the mean squared error, and the closer to 0, the closer the
reconstruction is to the original image. We can see that the difference in the MSE between the
two reduced versions of the images is different but, in the reconstruction, this is barely
noticeable. More notably however, is the ratio. This tells us how many times less data needs to
be stored for each image. This is significant due to the fact that we are reproducing near
identical images, however in the Rank (50) images we store approximately 15 times and 35
times less numbers respectively. This is even more significant for the Rank (20, 50, 10) images
where we store nearly 30 times and 65 times less than the original. 65 times less data while
retaining near identical image resolution is a great trade off.
Question 3a: Diode Current Modelling

To approximate the given piecewise function over the interval [-1,1], we first construct a
orthonormal basis of polynomials using the Gram-Schmidt process with respect to the inner
product
?
⟨𝑓, 𝑔⟩ = _ 𝑓(𝑥)𝑔(𝑥) 𝑑𝑥
@?

The basis consists of polynomials of increasing degree, each orthogonal to the others and
normalized to have unit length. Once the basis is constructed, the function 𝑓(𝑥) is projected
onto the basis to find the least squares polynomials approximation.
For a degree 3 approximation, the projection includes the constant, linear, quadratic and cubic
basis functions. This gives the polynomial that best approximates 𝑓(𝑥) in the least squares
sense up to degree 3.
For a degree 4 approximation, the quartic term is included as well, this improves the accuracy
by allowing the polynomial to better adapt to the curvature of 𝑠𝑖𝑛(𝜋𝑥) on [0,1].

When plotted alongside the original function, the degree 3 polynomial captures the general
sinusoidal shape on [0,1], but it less accurate near the transition at 𝑥 = 0. It slightly
underestimates the peak near 𝑥 = 1 and smooths out the sharp change in slope at the
transition. The degree 4 polynomial has a noticeably better fit. It follows the sinusoidal curve
more closely and reduced the deviation near the jump. Some small oscillation still occurs for
𝑥 < 0. This behavior is expected with this degree polynomial approximation.
The figure below shows the original function and both polynomial approximations plotted on
the same set of axes.

Figure 12: Degree 3 and 4 Polynomial Function Approximates of the Diode Function over the interval [-1, 1]
Question 3b:

The accuracy of the degree 4 polynomial approximation can be evaluated by plotting the
absolute error between the original function 𝑓(𝑥) and its approximation 𝑓2196 (𝑥). The error
function is defined as:
Error(𝑥) = d𝑓(𝑥) − 𝑓deg4 (𝑥)d.
When this error is plotted over the interval [-1,1], the largest error occurs at 𝑥 = 1, with a
notable error also visible at 𝑥 = 0. This error at 𝑥 = 0 arises because the original function
transition from a constant value (𝑓(𝑥) = 0 𝑓𝑜𝑟 𝑥 < 0) to a sinusoidal shape (𝑓(𝑥) =
𝑠𝑖𝑛(𝜋𝑥) for 𝑥 ≥ 0). Although the function is continuous at this point, the slope changes
abruptly. Polynomials are smooth functions and cannot perfectly reproduce sudden changes
in slope, leading to a noticeable error at 𝑥 = 0. However, the maximum error is found at the
end point 𝑥 = 1, where the polynomial approximation slightly overshoots the function - a
common error near interval boundaries.

Figure 13: Absolute error of a degree 4 polynomial approximation

Question 3c:

To achieve a more accurate approximation of the function 𝑓(𝑥), we want the minimum
polynomial degree required to reduce the maximum absolute error below 0.1 over the interval
[-1,1]. Using the orthonormal polynomial basis generated via the Gram-Schmidt process, we
computed successive polynomial approximations of increasing degree and evaluated their
maximum absolute errors.
The process can be implemented in SymPy with a simple condition:
target_error = 0.1
sample_points = [-1 + i*0.005 for i in range(401)]

for k in range(1, 21):

approx = sum([Link](func*v, (x, -1, 1)) * v for v
in orthonormal_basis[:k])

max_err = max(abs([Link](x, p) - [Link](x, p)) for p

in sample_points)

if max_err < target_error:

best_degree = k

Listing 1. Calculating degree polynomial with max error less than 0.1

Using this approach, it was found that a degree 10 polynomial is sufficient to reduce the
maximum error below 0.1. The figure below shows the absolute error between the original
function and the degree 10 polynomial approximation. As you can see the max error (at 𝑥 =
0) is less than 0.1 (~0.09).

Figure 14: Absolute error for the degree 10 approximation.

The graph below shows the original function 𝑓(𝑥) plotted alongside the degree 10
polynomial approximation. As expected, the higher degree polynomial provided a much
closer match to the sinusoidal portion on [0,1] and significantly recued the error near the
transition at 𝑥 = 0. The approximation follows the overall shape of the function closely
across the entire interval, with only small deviation remaining, consistent with the observed
maximum error being less than 0.1.

Figure 15: Comparison of the original function and the degree 10 polynomial approximation

Question 3d:

When approximating the diode current function, SymPy treats the function symbolically and
allows the construction of orthonormal polynomials using exact symbolic integration. By
using SymPy to perform the Gram–Schmidt process, the polynomials are generated
automatically without the need for manual integration or normalization, which can be very
time consuming for higher-degree polynomials. This approach ensures that the projection
coefficients are mathematically exact and the approximations retain the correct form across
the interval. SymPy also provides built-in plotting functionality, allowing both the original
function and the approximations to be visualised directly within the same environment.
However, symbolic computation can be slow, particularly as the degree of the polynomials
increases, and the resulting expressions can become complex and difficult to interpret. In
addition, SymPy is not optimized for numerical performance, so evaluating or plotting dense
symbolic expressions can be inefficient. In some cases, certain integrals may not simplify
automatically and require manual adjustment. Despite these limitations, SymPy provides an
exact and educational approach to polynomial approximation, making it ideal for
understanding the underlying mathematics even if it is less practical for large-scale or high-
degree computations.