forked from GiggleLiu/ModernScientificComputing
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path4.linearequation.jl
2082 lines (1738 loc) · 67.9 KB
/
4.linearequation.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
### A Pluto.jl notebook ###
# v0.19.22
using Markdown
using InteractiveUtils
# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
quote
local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
local el = $(esc(element))
global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
el
end
end
# ╔═╡ 37dc4b11-8849-41ad-9afe-5f792ea1fab8
using PlutoUI
# ╔═╡ 8dc256a3-06be-4df9-aa88-1ef1e574055e
using Test, LinearAlgebra
# ╔═╡ c5aaa7ec-9af4-45eb-86f1-e87cef9c8384
using BenchmarkTools
# ╔═╡ d959c85c-0895-475c-b9f8-57368fc5e744
using Plots
# ╔═╡ cdaafdfa-9252-410f-af0c-e47a145fbfe3
TableOfContents(; depth=2)
# ╔═╡ 2464e2fc-1a00-4572-b54c-f87b3239f9ad
html"<button onclick='present();'>present</button>"
# ╔═╡ f6fa6216-0945-4c98-b507-7c91600d8c2f
md"## About assignments"
# ╔═╡ a0fa8273-cf22-4cb4-ad4b-b83aee874556
md"1. You should submit your homework by opening a pull request to [our Github repo](https://github.com/GiggleLiu/ModernScientificComputing). Xuanzhao Gao will comment on your PR, then you should resolve the issues to get your PR merged. You should submit different PRs for different assignments.
2. We will grade based on your merged PRs. You will not be failed if your submition is complete.
4. You may check the answers of other students, but you must credit him in your PR description, e.g.
```
### Reference:
* PR #3
* PR #4
```
"
# ╔═╡ 78b6f5be-0aec-4e40-8512-fb58c253a7e9
md"# System of Linear Equations
Let $A\in \mathbb{R}^{n\times n}$ be a invertible square matrix and $b \in \mathbb{R}^n$ be a vector. Solving a linear equation means finding a vector $x\in\mathbb{R}^n$ such that
```math
A x = b
```
"
# ╔═╡ c1cf7bba-3d47-4c8c-8f1b-64808c300f72
md"Quiz: In a cage, chickens and rabbits add up to 35 heads and 94 feet. Please count the number of chickens and rabbits."
# ╔═╡ 24e421db-547b-416b-8c37-582e9aa78165
md"### Table of contents"
# ╔═╡ e5c4aeef-5818-406c-be55-b6b010c960db
md"* Gaussian elimination algorithm
* Pivoting technique
* Sensitivity analysis and condition number
* Getting matrix inverse with Gauss-Jordan elimination
* Solving linear equations for special matrices (optional)
"
# ╔═╡ a7203d10-cffd-4f08-80c2-f6b1465210bd
md"## Schetch of solving the linear equation"
# ╔═╡ 3a411c5f-264b-4f81-98ba-53eecdad3d28
md"""
One can solve a linear equation by following these steps:
1. Decompose the matrix $A \in \mathbb{R}^{n\times n}$ into $L \in \mathbb{R}^{n\times n}$ and $U \in \mathbb{R}^{n\times n}$ matrices using a method such as [Gaussian elimination](https://en.wikipedia.org/wiki/Gaussian_elimination) or [Crout's method](https://en.wikipedia.org/wiki/Crout_matrix_decomposition).
2. Rewrite the equation $Ax = b$ as $LUx = b$.
3. Solve for y in $Ly = b$ by [forward substitution](https://en.wikipedia.org/wiki/Triangular_matrix). This involves substituting the values of $y$ into the equation one at a time, starting with the first row and working downwards.
4. Solve for $x$ in $Ux = y$ by [back substitution](https://en.wikipedia.org/wiki/Triangular_matrix). This involves substituting the values of $x$ into the equation one at a time, starting with the last row and working upwards.
"""
# ╔═╡ 4015a2f9-a281-44ab-974c-7cd46041d105
md"""
## Solving tridiagonal linear equation
```math
Lx = b
```
"""
# ╔═╡ 53e51d33-9d03-4473-8ed8-aba1ef9ef105
md"""
where $L \in \mathbb{R}^{n\times n}$ is a lower triangular matrix defined as
```math
L = \left(\begin{matrix}
l_{11} & 0 & \ldots & 0\\
l_{21} & l_{22} & \ldots & 0\\
\vdots & \vdots & \ddots & \vdots\\
l_{n1} & l_{n2} & \ldots & l_{nn}
\end{matrix}\right)
```
"""
# ╔═╡ 02605f8e-b457-11ed-19c4-7d6d14b67507
md"""## Algorithm: Forward-Substitution for Lower Triangular System"""
# ╔═╡ c86b3ac8-3ae2-40bb-8988-d9b9250d627d
md"""
Forward substitution is an algorithm used to solve a system of linear equations with a lower triangular matrix. It involves solving for the unknowns in a forward direction, starting from the first equation and moving towards the last. Here's an example of forward substitution algorithm in the matrix form:
Consider the following system of lower triangular linear equations:
```math
L = \left(\begin{matrix}
3 & 0 & 0\\
2 & 5 & 0\\
1 & 4 & 2
\end{matrix}\right)
\left(\begin{matrix}
x_1\\
x_2\\
x_3
\end{matrix}\right) =
\left(\begin{matrix}
9\\
12\\
13
\end{matrix}\right)
```
To solve for $x_1$, $x_2$, and $x_3$ using forward substitution, we start with the first equation:
```math
3x_1 + 0x_2 + 0x_3 = 9
```
Solving for $x_1$, we get:
```math
x_1 = 3
```
Substituting $x = 3$ into the second equation (row), we get:
```math
2(3) + 5x_2 + 0x_3 = 12
```
Solving for $x_2$, we get:
```math
x_2 = (12 - 6) / 5 = 1.2
```
Substituting $x = 3$ and $x_2 = 1.2$ into the third equation (row), we get:
```math
1(3) + 4(1.2) + 2x_3 = 13
```
Solving for $x_3$, we get:
```math
x_3 = (13 - 3 - 4(1.2)) / 2 = 1.5
```
Therefore, the solution to the system of equations is:
```math
x = \left(\begin{matrix}\
3\\
1.2\\
1.5
\end{matrix}\right)
```
"""
# ╔═╡ fb583209-eb00-41ce-9fc6-906898327ecd
md"""
It can be summarized to the following algorithm
```math
x_1 = b_1/l_{11},~~~ x_i = \left(b_i - \sum_{j=1}^{i-1}l_{ij}x_j\right)/l_{ii},~~ i=2, ..., n
```
"""
# ╔═╡ 35d1d353-c3bd-4e2d-90e0-e3031a2133c2
md"### The implementation"
# ╔═╡ 5c73b504-f8f2-429a-99a2-5f5d98f7e1e6
md"""
We implement the above algorithm in Julia language.
"""
# ╔═╡ 742a79b6-6ced-467e-961b-814da8b65af7
function back_substitution!(l::AbstractMatrix, b::AbstractVector)
n = length(b)
@assert size(l) == (n, n) "size mismatch"
x = zero(b)
# loop over columns
for j = 1:n
# stop if matrix is singular
if iszero(l[j, j])
error("The lower triangular matrix is singular!")
end
# compute solution component
x[j] = b[j] / l[j, j]
for i = j+1:n
# update right hand side
b[i] = b[i] - l[i, j] * x[j]
end
end
return x
end
# ╔═╡ e5e54678-4d3a-40ae-8b39-d5c50a0617b2
md"""
We can write a test for this algorithm.
"""
# ╔═╡ d3e1ef10-3f99-4561-87db-47c0fc47e9e9
@testset "back substitution" begin
# create a random lower triangular matrix
l = tril(randn(4, 4))
# target vector
b = randn(4)
# solve the linear equation with our algorithm
x = back_substitution!(l, copy(b))
@test l * x ≈ b
# The Julia's standard library `LinearAlgebra` contains a native implementation.
x_native = LowerTriangular(l) \ b
@test l * x_native ≈ b
end
# ╔═╡ 7649913c-523f-4fa3-bf1e-ea1566b55114
md"## The LU decomposition with Gaussian elimination
LU decomposition is a method for solving linear equations that involves breaking down a matrix into lower and upper triangular matrices. The $LU$ decomposition of a matrix $A$ is represented as $A = LU$, where $L$ is a lower triangular matrix and $U$ is an upper triangular matrix.
"
# ╔═╡ 24fa07b8-6a9c-47d6-8f72-9346f1325ab4
md"## The elementary elimination matrix"
# ╔═╡ 851510d9-7703-4545-8c00-ec94d8735a85
md"""
An elementary elimination matrix is a matrix that is used in the process of Gaussian elimination to transform a system of linear equations into an equivalent system that is easier to solve. It is a square matrix that is obtained by performing a single elementary row operation on the identity matrix.
"""
# ╔═╡ af1f932e-a6db-41ce-ace3-170228bb183b
# ```math
# (M_k)_{ij} = \begin{cases}
# \delta_{ij} & i= j,\\
# - a_{ik}/a_{kk} & i > j \land j = k, \\
# 0 & {\rm otherwise}.
#\end{cases}
#```
md"""
Let $A = (a_{ij})$ be a square matrix of size $n \times n$. The $k$th elementary elimination matrix for it is defined as
```math
M_k = \left(\begin{matrix}
1 & \ldots & 0 & 0 & 0 & \ldots & 0\\
\vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots\\
0 & \ldots & 1 & 0 & 0 & \ldots & 0\\
0 & \ldots & 0 & 1 & 0 & \ldots & 0\\
0 & \ldots & 0 & -m_{k+1} & 1 & \ldots & 0\\
\vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots\\
0 & \ldots & 0 & -m_{n} & 0 & \ldots & 1\\
\end{matrix}\right)
```
where $m_i = a_{ik}/a_{kk}$.
"""
# ╔═╡ a890cdb6-34b6-4a3c-8749-e1543e9764ba
md"""
By applying this elimentary elimination matrix $M_1$ on $A$, we can obtain a new matrix with the $a_{i1}' = 0$ for all $i>1$.
```math
M_1 A = \left(\begin{matrix}
a_{11} & a_{12} & a_{13} & \ldots & a_{1n}\\
0 & a_{22}' & a_{23}' & \ldots & a_{2n}'\\
0 & a_{32}' & a_{33}' & \ldots & a_{3n}'\\
\vdots & \vdots & \vdots & \ddots & \vdots\\
0 & a_{n2}' & a_{n3}' & \ldots & a_{nn}'\\
\end{matrix}\right)
```
"""
# ╔═╡ 9640cc8d-f8c0-4d3b-b51a-3ad3865630bc
md"For $k=1,2,\ldots,n$, apply $M_k$ on $A$. We will have an upper triangular matrix.
```math
U = M_{n-1}\ldots M_1 A
```
Since $M_k$ is reversible (as we will show below), we have
```math
\begin{align}
&A = LU\\
&L = M_1^{-1} M_2^{-1} \ldots M_{n-1}^{-1},
\end{align}
```
"
# ╔═╡ 83e0cdc0-70d3-48d2-bf1c-8ed364652aed
md"## Properties of elementary elimination matrices
1. Its inverse can be computed in $O(n)$ time
```math
M_k^{-1} = 2I - M_k
```
2. The multiplication of two elementary matrices can be computed in $O(n)$ time
```math
M_k M_{k' > k} = M_k + M_{k'} - I
```
"
# ╔═╡ 2b3596d7-17a3-46bb-b687-302c1255a4af
md"## Coding: Properties of elimination matrices"
# ╔═╡ b8cf4822-d91b-47f4-86e4-a36d668b5256
A3 = [1 2 2; 4 4 2; 4 6 4]
# ╔═╡ faa03439-b526-4c27-a9dc-593b5a32cb04
function elementary_elimination_matrix(A::AbstractMatrix{T}, k::Int) where T
n = size(A, 1)
@assert size(A, 2) == n
# create Elementary Elimination Matrices
M = Matrix{Float64}(I, n, n)
for i=k+1:n
M[i, k] = -A[i, k] ./ A[k, k]
end
return M
end
# ╔═╡ 416f9a74-2d88-4788-92b1-847eda6539cc
md"The elimination"
# ╔═╡ b8379e6d-768d-4a43-8370-1583324750f3
elementary_elimination_matrix(A3, 1)
# ╔═╡ e478efbc-a75f-4fb4-8ea8-f24ea7f8a9c2
elementary_elimination_matrix(A3, 1) * A3
# ╔═╡ 73fa4829-33dc-46e7-935d-b4c78bdb02d3
md"The inverse"
# ╔═╡ fc286cee-c912-4395-8af6-77702d5fef31
inv(elementary_elimination_matrix(A3, 1))
# ╔═╡ b4282280-a825-456f-90ec-cfc274c68864
md"The multiplication"
# ╔═╡ 32ce956f-a90f-44fb-8473-438dd2a6ba78
elementary_elimination_matrix(A3, 2)
# ╔═╡ 2a8da2fe-f21e-44ac-9253-583c06dd5134
inv(elementary_elimination_matrix(A3, 1)) * inv(elementary_elimination_matrix(A3, 2))
# ╔═╡ 87de48b5-4d4e-4b8a-9f2c-4b6dc328a502
md"## LU Factorization by Gaussian Elimination"
# ╔═╡ b27b9f35-bac1-4d4c-9543-3733ed21e039
md"""
A naive implementation of elimentary elimination matrix is as follows
"""
# ╔═╡ a4b32a01-2d7f-4933-a3cb-2ab2d0c6cb2e
function lufact_naive!(A::AbstractMatrix{T}) where T
n = size(A, 1)
@assert size(A, 2) == n
M = Matrix{T}(I, n, n)
for k=1:n-1
m = elementary_elimination_matrix(A, k)
M = M * inv(m)
A .= m * A
end
return M, A
end
# ╔═╡ 4d8276dc-9356-47f8-afa6-5329dde63468
lufact_naive!(copy(A3))
# ╔═╡ 405450f7-807d-45a8-94ac-54ef9972b234
@testset "naive LU factorization" begin
A = [1 2 2; 4 4 2; 4 6 4]
L, U = lufact_naive!(copy(A))
@test L * U ≈ A
end
# ╔═╡ 1df2f386-2b55-4fb0-a035-64cf454e2698
md"The above implementation has time complexity $O(n^4)$ since we did not use the sparsity of elimentary elimination matrix. A better implementation that gives $O(n^3)$ time complexity is as follows."
# ╔═╡ fadba4ba-efca-48f7-bf31-959f458f64d5
function lufact!(a::AbstractMatrix)
n = size(a, 1)
@assert size(a, 2) == n "size mismatch"
m = zero(a)
m[1:n+1:end] .+= 1
# loop over columns
for k=1:n-1
# stop if pivot is zero
if iszero(a[k, k])
error("Gaussian elimination fails!")
end
# compute multipliers for current column
for i=k+1:n
m[i, k] = a[i, k] / a[k, k]
end
# apply transformation to remaining sub-matrix
for j=k+1:n
for i=k+1:n
a[i,j] -= m[i,k] * a[k, j]
end
end
end
return m, triu!(a)
end
# ╔═╡ f2bab70d-847f-495b-ba30-409156a9a874
lufact(a::AbstractMatrix) = lufact!(copy(a))
# ╔═╡ 480b6065-f9ca-4770-baa4-f1451fc046b8
@testset "LU factorization" begin
a = randn(4, 4)
L, U = lufact(a)
@test istril(L)
@test istriu(U)
@test L * U ≈ a
end
# ╔═╡ b81714d1-9a3c-4800-b478-68cabbc5d10d
A4 = randn(4, 4)
# ╔═╡ 1f1267d0-9ff2-44a4-b418-a42953159a1e
lufact(A4)
# ╔═╡ 5e426a5f-65cf-4adf-8b2b-ea28534e529a
md"Julia language has a much better implementation in the standard library `LinearAlgebra`."
# ╔═╡ 9d30b63e-f437-420b-8c22-037d34e1be2f
julia_lures = lu(A4, NoPivot()) # the version we implemented above has no pivot
# ╔═╡ 4a18900a-9051-4451-9736-724b19c3bd5d
julia_lures.U
# ╔═╡ aacd052d-b6e7-4f40-b08a-872aa9420c46
typeof(julia_lures)
# ╔═╡ c86a4d3c-adf3-429a-a4b4-8a0aaec83b43
fieldnames(julia_lures |> typeof)
# ╔═╡ 2fda47ad-08f2-45c5-8f47-0eec1aa49aac
md"## Issue: how to handle small diagonal entries?"
# ╔═╡ 0d8dc065-4487-41bd-a80d-45f1fc1bf550
md"""
The above Gaussian elimination process is not stable if any diagonal entry in $A$ has a value that close to zero.
"""
# ╔═╡ 77396c1b-f23d-41da-b9e2-d6980a090217
@bind ϵ Slider(-20:0.01:0.0; default=-2, show_value=true)
# ╔═╡ 97e373df-32cc-491e-8066-e3a33b9b2764
small_diagonal_matrix = [10^(ϵ) 1; 1 1]
# ╔═╡ 1201b122-b071-4f13-83e0-6bec0de08e2b
lures = lufact(small_diagonal_matrix)
# ╔═╡ 7a7e8ac0-2af9-42ec-ab17-54d4f86baf70
md"This issue is not always related to the rank deficiency (or ill-condition). A remedy to this issue is to permute the rows of $A$ before factorizing it.
For example:"
# ╔═╡ a4c024e4-83a7-4356-9b8e-168ce9ee1cbf
lufact(small_diagonal_matrix[end:-1:1, :])
# ╔═╡ 99381db6-0735-4bfe-9d20-1e16e08f13c3
md"This technique is called pivoting."
# ╔═╡ 54db9443-a0f9-48f1-97bb-bcab2ae3cfd5
md"""
## Pivoting technique
"""
# ╔═╡ 178c4193-a45d-460f-a22e-b15f6604b5dd
md"""
LU factoriazation (or Gaussian elimination) with pivoting is defined as
```math
P A = L U
```
where $P$ is a permutation matrix.
"""
# ╔═╡ 5a034215-a5da-4fea-8dc3-004310f5e7be
md"Pivoting in Gaussian elimination is the process of selecting a pivot element in a matrix and then using it to eliminate other elements in the same column or row. The pivot element is chosen as the largest absolute value in the column, and its row is swapped with the row containing the current element being eliminated if necessary. This is done to avoid division by zero or numerical instability, and to ensure that the elimination process proceeds smoothly. Pivoting is an important step in Gaussian elimination, as it ensures that the resulting matrix is in reduced row echelon form and that the solution to the system of equations is accurate."
# ╔═╡ 9bd48ccc-549a-4057-aa91-16ad79536583
md"""
## Gaussian Elimination process with Partial Pivoting
Let $A=(a_{ij})$ be a square matrix of size $n\times n$. The Gaussian elimination process with partial pivoting can be represented as
```math
M_{n-1}P_{n-1}\ldots M_2P_2M_1P_1 A = U
```
"""
# ╔═╡ 95c78df3-2283-4230-8d8d-9996e4d6fe36
md"Here we emphsis that $P_{k}$ and $M_{j<k}$ commute."
# ╔═╡ e9ba7dc4-1494-4601-8eb7-72127e92815d
md"## LU Factoriazation by Gaussian Elimination with Partial Pivoting"
# ╔═╡ 5a8e40b0-1cd7-4a71-94da-64cfd9ab9755
md"""
A Julia implementation of the Gaussian elimination with partial pivoting is
"""
# ╔═╡ 9b4fa87d-c40c-4677-82f2-2fff3b87f4c6
function lufact_pivot!(a::AbstractMatrix)
n = size(a, 1)
@assert size(a, 2) == n "size mismatch"
m = zero(a)
P = collect(1:n)
# loop over columns
@inbounds for k=1:n-1
# search for pivot in current column
val, p = findmax(x->abs(a[x, k]), k:n)
p += k-1
# find index p such that |a_{pk}| ≥ |a_{ik}| for k ≤ i ≤ n
if p != k
# swap rows k and p of matrix A
for col = 1:n
a[k, col], a[p, col] = a[p, col], a[k, col]
end
# swap rows k and p of matrix M
for col = 1:k-1
m[k, col], m[p, col] = m[p, col], m[k, col]
end
P[k], P[p] = P[p], P[k]
end
if iszero(a[k, k])
# skip current column if it's already zero
continue
end
# compute multipliers for current column
m[k, k] = 1
for i=k+1:n
m[i, k] = a[i, k] / a[k, k]
end
# apply transformation to remaining sub-matrix
for j=k+1:n
akj = a[k, j]
for i=k+1:n
a[i,j] -= m[i,k] * akj
end
end
end
m[n, n] = 1
return m, triu!(a), P
end
# ╔═╡ 13493579-9189-4c9a-9361-430868427b59
@testset "lufact with pivot" begin
n = 5
A = randn(n, n)
L, U, P = lufact_pivot!(copy(A))
pmat = zeros(Int, n, n)
setindex!.(Ref(pmat), 1, 1:n, P)
@test L ≈ lu(A).L
@test U ≈ lu(A).U
@test pmat * A ≈ L * U
end
# ╔═╡ f3a506b5-669f-4e8f-9424-e04532b9a87a
md"""
If you are interested in knowing the performance of our implementation, please check the following check box.
"""
# ╔═╡ d1b7f00e-2b7e-4c15-8231-a816c59472f9
@bind benchmark_lu CheckBox()
# ╔═╡ 3bcf8d78-3678-4749-99ac-0a7e43830d2d
if benchmark_lu let
n = 200
A = randn(n, n)
@benchmark lufact_pivot!($A)
end end
# ╔═╡ e206e554-64cd-422d-963d-8284eb9625fe
if benchmark_lu let
n = 200
A = randn(n, n)
@benchmark lu($A)
end end
# ╔═╡ ed52cf76-ed66-44af-b98b-b58551bd217f
md"""
## Complete pivoting
The complete pivoting also allows permuting columns. It produces better numerical stability but is also harder to implement. In most practical using cases, partial pivoting is good enough.
```math
P A Q = L U
```
"""
# ╔═╡ e113247a-b2a5-4acc-8714-342ecb191ac4
md"# Sensitivity Analysis"
# ╔═╡ 9231de22-d474-4548-b984-fcc3dda27113
md"""
Sensitivity analysis in linear algebra is the study of how changes in the input data or parameters of a linear system affect the output or solution of the system. It involves analyzing the sensitivity of the solution to changes in the coefficients of the system, the right-hand side vector, or the constraints of the problem. Sensitivity analysis can be used to determine the effect of small changes in the input data on the optimal solution, to identify critical parameters that affect the solution, and to assess the robustness of the solution to uncertainties or variations in the input data.
"""
# ╔═╡ 25f870c1-e70a-48b1-a0ce-4c2b9a523d60
md"## Issue: An Ill Conditioned Matrix"
# ╔═╡ 0c5a1e3f-7383-4705-b63e-e64f4a513f85
md"""
An ill conditioned matrix may produce unreliable result, or the output is very sensitive to the input. The following is an example of a matrix close to singular.
"""
# ╔═╡ b333fc58-0981-4d0f-ba49-0ee15cc78aad
md"""
```math
A = \left(\begin{matrix}
0.913 & 0.659\\
0.457 & 0.330
\end{matrix}
\right)
```
"""
# ╔═╡ 0beab22f-5411-403c-a116-cf552083382e
ill_conditioned_matrix = [0.913 0.659; 0.457 0.330]
# ╔═╡ 63cddfd1-bdef-4601-ab57-48314952fb73
lures2 = lufact(ill_conditioned_matrix)
# ╔═╡ c5b109da-ab9c-4a44-ab2d-168b9711c2cc
lures2
# ╔═╡ 713267d1-08f0-4c41-9523-46128ce40430
cond(ill_conditioned_matrix)
# ╔═╡ 51433b0f-8fce-4241-bb84-bcea9a687abe
md"## The relative error"
# ╔═╡ 26207f81-1f6e-4708-8c04-422a1133c5ac
md"""
The relevant error in floating number system is the relative error.
* Absolute error: $\|x - \hat x\|$
* Relative error: $\frac{\|x - \hat x\|}{\|x\|}$
where $\|\cdot\|$ is a measure of size.
"""
# ╔═╡ 025101f0-a1c4-47da-bff0-c0ae0cd0da6e
md"## Numeric experiment: Floating point numbers have \"constant\" relative error"
# ╔═╡ 424638de-6651-4fe2-bc44-82f9675dc43d
eps(Float64)
# ╔═╡ 26978647-a790-4ba2-be21-60251c47d919
let
n = 1000
reltol = zeros(2n+1)
for i=-n:n
f = 2.0^i
reltol[i+n+1] = log10(eps(f)) - log10(f)
end
plot(-n:n, reltol; label="relative error")
end
# ╔═╡ c10acde5-e22a-4f79-be4f-f7b3bbb0ece9
eps(1.0)/1.0
# ╔═╡ 8e055682-baef-4c78-b044-d22d51676844
eps(2.0)/2.0
# ╔═╡ ee1aba8e-5992-4ff1-9ae1-09d35a2abf00
eps(sqrt(2))/sqrt(2)
# ╔═╡ d1c29e2e-e67d-41d3-bbff-9f384cf0c84d
md"## (Relative) Condition Number"
# ╔═╡ 2f789e07-69a9-4786-8d82-e363afa4323f
md"
Condition number is a measure of the sensitivity of a mathematical problem to changes or errors in the input data. It is a way to quantify how much the output of a function or algorithm can vary due to small changes in the input. A high condition number indicates that the problem is ill-conditioned, meaning that small errors in the input can lead to large errors in the output. A low condition number indicates that the problem is well-conditioned, meaning that small errors in the input have little effect on the output.
In short, the (relative) condition number of an operation $f$ with input $x$ measures the relative error amplication power of $f$ with input $x$, which is formally defined as
"
# ╔═╡ df17edcb-8b98-4c2e-a628-1a1381b7cf3e
md"""
```math
\lim _{\varepsilon \rightarrow 0^{+}}\,\sup _{\|\delta x\|\,\leq \,\varepsilon }{\frac {\|\delta f(x)\|/\|f(x)\|}{\|\delta x\|/\|x\|}}.
```
"""
# ╔═╡ 46176c75-f23e-4aef-a54b-58aa43941855
md"Quiz 1: What is the condition number of"
# ╔═╡ fe41dd53-f69f-4382-bde1-af86a3a92aef
md"""
```math
y = \exp(x)
```
"""
# ╔═╡ 8651c71c-eb8b-467b-9c9c-f2a173e5b91d
md"Quiz 2: What is the condition number of"
# ╔═╡ d9d88d11-0dd0-4e80-b476-e9b008dfb84c
md"""
```math
a + b
```
"""
# ╔═╡ 88ad6fff-04d8-477d-a664-c1a6522964cc
md"""
With the obtained result, show why we should avoid subtracting two big floating point numbers?
"""
# ╔═╡ bb5e9e72-8cc7-4c0c-96c2-9c1ca9d1b38c
md"## Measuring the size vectors and matrices"
# ╔═╡ e1026249-5d13-418b-a908-b6dc77175434
md"The vector $p$-norm is formally defined as"
# ╔═╡ f256bf0d-e2bb-49f4-882f-12508e68ff07
md"""
```math
\|v\|_p = \left(\sum_i{|v_i|^p}\right)^{1/p}
```
"""
# ╔═╡ f2da27a2-af5f-4f2a-a6d5-fb98b8dcd261
md"""Similarly, the matrix $p$-norm is formally defined as
```math
\|A\|_p = \max_{x\neq 0} \frac{\|Ax\|_p}{\|x\|_p}
```
"""
# ╔═╡ 49d62ec4-bce0-4620-93cc-278f149f848f
md"## Examples: Vector and Matrix Norms"
# ╔═╡ 33a6261a-2b21-4b37-bcd8-b5c0db6b509f
norm([3, 4], 2)
# ╔═╡ 1956af2c-c88f-483d-8e1c-8ed4004e996d
norm([3, 4], 1)
# ╔═╡ 61c1c35c-f7be-495b-90f5-9a3772f254f5
norm([3, 4], Inf)
# ╔═╡ 43731510-dfe7-48c8-8a4d-d8313d1d762b
# l0 norm is not a true norm
norm([3, 4], 0)
# ╔═╡ a0af42b2-1bba-423c-b5f7-0a9d1a2adbb1
norm([3, 0], 0)
# ╔═╡ 79bc3efb-61a1-4241-b787-2726fb2a4a97
mat = randn(2, 2)
# ╔═╡ 9854ad79-d06e-4d60-86d5-54ade042d49a
opnorm(mat, 1)
# ╔═╡ 4628cb4b-649f-4e84-a6d5-e50109f8cc6c
opnorm(mat, Inf)
# ╔═╡ 424f63ce-f5b2-456f-8c4a-63f18877d730
opnorm(mat, 2)
# ╔═╡ 71e6ed15-bb0f-4510-901c-c9f53b3e701c
opnorm(mat, 0)
# ╔═╡ 898a511d-7c18-4eb0-b68b-3c1dd8208cfe
cond(mat)
# ╔═╡ 5653168f-80b5-4435-a711-d8a1d8f9c429
md"## Condition Number of a Linear Operator"
# ╔═╡ c165a6b1-78da-40b5-b925-a657193ade6a
md"""
The condition number of a linear system $Ax = b$ is defined as
"""
# ╔═╡ a1a55972-4baa-4a29-8c3f-ac7e13036136
md"""
```math
{\rm cond}(A) = \|A\| \|A^{-1}\|
```
"""
# ╔═╡ 0b1c4d61-b748-4e9d-9a29-9ee50fa90e79
md"""
Using the defintion of condition number, we have
```math
\begin{align}
{\rm cond(A, x)}&=\lim _{\varepsilon \rightarrow 0^{+}}\,\sup _{\|\delta x\|\,\leq \,\varepsilon }{\frac {\|\delta (Ax)\|/\|A x\|}{\|\delta x\|/\|x\|}}\\
&=\lim _{\varepsilon \rightarrow 0^{+}}\,\sup _{\|\delta x\|\,\leq \,\varepsilon }{\frac {\|A\delta x\|\|x\|}{\|\delta x\|\|Ax\|}}
\end{align}
```
Let $y = Ax$, we have
```math
\begin{align}
{\rm cond(A, x)}&=\lim _{\varepsilon \rightarrow 0^{+}}\,\sup _{\|\delta x\|\,\leq \,\varepsilon }{\frac {\|A\delta x\|\|A^{-1}y\|}{\|\delta x\|\|y\|}}\\
&=\|A\|\frac{\|A^{-1}y\|}{\|y\|}
\end{align}
```
Suppose we want to get an upper bound for any input $x$, then using the definiton of matrix norm, we have
```math
{\rm cond}(A) = \|A\|\sup_y \frac{\|A^{-1}y\|}{\|y\|} = \|A\| \|A^{-1}\|
```
"""
# ╔═╡ 8ee2c9f5-3dbf-4ab0-a638-5a8c2d5d1367
md"""## Numeric experiment: Numeric experiment on condition number"""
# ╔═╡ 180b4de9-9b37-47c2-905a-0156210be5cc
md"""
We randomly generate matrices of size $10\times 10$ and show the condition number approximately upper bounds the numeric error amplification factor of a linear equation solver.
"""
# ╔═╡ b71889e3-85d0-4873-a94a-63a9c4c0c009
let
n = 1000
p = 2
errors = zeros(n)
conds = zeros(n)
for k = 1:n
A = rand(10, 10)
b = rand(10)
dx = A \ b
sx = Float32.(A) \ Float32.(b)
errors[k] = (norm(sx - dx, p)/norm(dx, p)) / (norm(b-Float32.(b), p)/norm(b, p))
conds[k] = cond(A, p)
end
plt = plot(conds, conds; label="condition number", xlim=(1, 10000), ylim=(1, 10000), xscale=:log10, yscale=:log10)
scatter!(plt, conds, errors; label="samples")
end
# ╔═╡ 927272dc-d60c-4263-8139-7f43887df52c
md"## Positive definite symmetric matrix"
# ╔═╡ eac7780a-5e0b-4518-907f-9ff75dfb4b2f
md"""
* (Real) Symmetric: $A = A^T$,
* Positive definite: $x^T A x > 0$ for all $x \neq 0$.
"""
# ╔═╡ 8dced7f2-8eab-41d8-95f7-5c23e0035466
md"## Cholesky decomposition
Cholesky decomposition is a method of decomposing a positive-definite matrix into a product of a lower triangular matrix and its transpose. It is named after the mathematician André-Louis Cholesky, who developed the method in the early 1900s. The Cholesky decomposition is useful in many areas of mathematics and science, including linear algebra, numerical analysis, and statistics. It is often used in solving systems of linear equations, computing the inverse of a matrix, and generating random numbers with a given covariance matrix. The Cholesky decomposition is computationally efficient and numerically stable, making it a popular choice in many applications.
Given a positive definite symmetric matrix $A\in \mathbb{R}^{n\times n}$, the Cholesky decomposition is formally defined as
```math
A = LL^T,
```
where $L$ is an upper triangular matrix.
"
# ╔═╡ d0bf8157-9139-418f-9bd8-2791e52aed8e
md"The implementation of Cholesky decomposition is similar to LU decomposition."
# ╔═╡ 79a5bc5d-a2bf-485a-95ee-07d3b8546e75
function chol!(a::AbstractMatrix)
n = size(a, 1)
@assert size(a, 2) == n
for k=1:n
a[k, k] = sqrt(a[k, k])
for i=k+1:n
a[i, k] = a[i, k] / a[k, k]
end
for j=k+1:n
for i=k+1:n
a[i,j] = a[i,j] - a[i, k] * a[j, k]
end
end
end
return a
end
# ╔═╡ 7945c565-d226-437c-83ad-5b1e3149c076
@testset "cholesky" begin
n = 10
Q, R = qr(randn(10, 10))
a = Q * Diagonal(rand(10)) * Q'
L = chol!(copy(a))
@test tril(L) * tril(L)' ≈ a
# cholesky(a) in Julia
end
# ╔═╡ 011e1b39-a73d-4eaf-8847-ae5fa5ffd2f0
md"""# Assignments
1. Get the relative condition number of division operation $a/b$.
2. Classify each of the following matrices as well-conditioned or ill-conditioned:
```math
(a). ~~\left(\begin{matrix}10^{10} & 0\\ 0 & 10^{-10}\end{matrix}\right)
```
```math
(b). ~~\left(\begin{matrix}10^{10} & 0\\ 0 & 10^{10}\end{matrix}\right)
```
```math
(c). ~~\left(\begin{matrix}10^{-10} & 0\\ 0 & 10^{-10}\end{matrix}\right)
```
```math
(d). ~~\left(\begin{matrix}1 & 2\\ 2 & 4\end{matrix}\right)
```
3. Implement the Gauss-Jordan elimination algorithm to compute matrix inverse. In the following example, we first create an augmented matrix $(A | I)$. Then we apply the Gauss-Jordan elimination matrices on the left. The final result is stored in the augmented matrix as $(I, A^{-1})$.

Task: Please implement a function `gauss_jordan` that computes the inverse for a matrix at any size. Please also include the following test in your submission.
```julia
@testset "Gauss Jordan" begin
n = 10
A = randn(n, n)
@test gauss_jordan(A) * A ≈ Matrix{Float64}(I, n, n)
end
```
"""
# ╔═╡ ca87e541-b520-4231-b949-4675456f6c0b
md"""Simular to computing Guassian elimination with elementary elimination matrices, computing the inverse can be done by repreatedly applying the Guass-Jordan elimination matrix and convert the target matrix to identity.
```math
SN_{n}N_{n-1}\ldots N_1 A = I
```
Then
```math
A^{-1} = SN_{n}N_{n-1}\ldots N_1
```
"""
# ╔═╡ c50dbe8c-6c09-47aa-a798-0673875e1e42
md"""
Here, the Gauss-Jordan elimination matrix $N_k$ eliminates the $k$th column except the diagonal element $a_{kk}$.
```math
N_k = \left(\begin{matrix}
1 & \ldots & 0 & -m_1 & 0 & \ldots & 0\\
\vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots\\
0 & \ldots & 1 & -m_{k-1} & 0 & \ldots & 0\\
0 & \ldots & 0 & 1 & 0 & \ldots & 0\\
0 & \ldots & 0 & -m_{k+1} & 1 & \ldots & 0\\
\vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots\\
0 & \ldots & 0 & -m_{n} & 0 & \ldots & 1\\
\end{matrix}\right)
```
where $m_i = a_i/a_k$.
"""
# ╔═╡ b429c34b-122e-4d3d-b8bd-96dab8b18368
md"""
$S$ is a diagonal matrix.
"""
# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
[compat]
BenchmarkTools = "~1.3.2"
Plots = "~1.38.5"
PlutoUI = "~0.7.50"
"""
# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised
julia_version = "1.8.5"
manifest_format = "2.0"
project_hash = "858d441359aea71f7f7b3a457de6e8128433ac4c"
[[deps.AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "8eaf9f1b4921132a4cff3f36a1d9ba923b14a481"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.4"
[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
version = "1.1.1"
[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
[[deps.BenchmarkTools]]
deps = ["JSON", "Logging", "Printf", "Profile", "Statistics", "UUIDs"]
git-tree-sha1 = "d9a9701b899b30332bbcb3e1679c41cce81fb0e8"
uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
version = "1.3.2"
[[deps.BitFlags]]
git-tree-sha1 = "43b1a4a8f797c1cddadf60499a8a077d4af2cd2d"
uuid = "d1d4a3ce-64b1-5f1a-9ba4-7e7e69966f35"
version = "0.1.7"