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
|
scipy.linalg
============
``scipy``\ ’s ``linalg`` module contains two functions,
``solve_triangular``, and ``cho_solve``. The functions can be called by
prepending them by ``scipy.linalg.``.
1. `scipy.linalg.solve_cho <#cho_solve>`__
2. `scipy.linalg.solve_triangular <#solve_triangular>`__
cho_solve
---------
``scipy``:
https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.cho_solve.html
Solve the linear equations
:raw-latex:`\begin{equation}
\mathbf{A}\cdot\mathbf{x} = \mathbf{b}
\end{equation}`
given the Cholesky factorization of :math:`\mathbf{A}`. As opposed to
``scipy``, the function simply takes the Cholesky-factorised matrix,
:math:`\mathbf{A}`, and :math:`\mathbf{b}` as inputs.
.. code::
# code to be run in micropython
from ulab import numpy as np
from ulab import scipy as spy
A = np.array([[3, 0, 0, 0], [2, 1, 0, 0], [1, 0, 1, 0], [1, 2, 1, 8]])
b = np.array([4, 2, 4, 2])
print(spy.linalg.cho_solve(A, b))
.. parsed-literal::
array([-0.01388888888888906, -0.6458333333333331, 2.677083333333333, -0.01041666666666667], dtype=float64)
solve_triangular
----------------
``scipy``:
https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve_triangular.html
Solve the linear equation
:raw-latex:`\begin{equation}
\mathbf{a}\cdot\mathbf{x} = \mathbf{b}
\end{equation}`
with the assumption that :math:`\mathbf{a}` is a triangular matrix. The
two position arguments are :math:`\mathbf{a}`, and :math:`\mathbf{b}`,
and the optional keyword argument is ``lower`` with a default value of
``False``. ``lower`` determines, whether data are taken from the lower,
or upper triangle of :math:`\mathbf{a}`.
Note that :math:`\mathbf{a}` itself does not have to be a triangular
matrix: if it is not, then the values are simply taken to be 0 in the
upper or lower triangle, as dictated by ``lower``. However,
:math:`\mathbf{a}\cdot\mathbf{x}` will yield :math:`\mathbf{b}` only,
when :math:`\mathbf{a}` is triangular. You should keep this in mind,
when trying to establish the validity of the solution by back
substitution.
.. code::
# code to be run in micropython
from ulab import numpy as np
from ulab import scipy as spy
a = np.array([[3, 0, 0, 0], [2, 1, 0, 0], [1, 0, 1, 0], [1, 2, 1, 8]])
b = np.array([4, 2, 4, 2])
print('a:\n')
print(a)
print('\nb: ', b)
x = spy.linalg.solve_triangular(a, b, lower=True)
print('='*20)
print('x: ', x)
print('\ndot(a, x): ', np.dot(a, x))
.. parsed-literal::
a:
array([[3.0, 0.0, 0.0, 0.0],
[2.0, 1.0, 0.0, 0.0],
[1.0, 0.0, 1.0, 0.0],
[1.0, 2.0, 1.0, 8.0]], dtype=float64)
b: array([4.0, 2.0, 4.0, 2.0], dtype=float64)
====================
x: array([1.333333333333333, -0.6666666666666665, 2.666666666666667, -0.08333333333333337], dtype=float64)
dot(a, x): array([4.0, 2.0, 4.0, 2.0], dtype=float64)
With get the same solution, :math:`\mathbf{x}`, with the following
matrix, but the dot product of :math:`\mathbf{a}`, and
:math:`\mathbf{x}` is no longer :math:`\mathbf{b}`:
.. code::
# code to be run in micropython
from ulab import numpy as np
from ulab import scipy as spy
a = np.array([[3, 2, 1, 0], [2, 1, 0, 1], [1, 0, 1, 4], [1, 2, 1, 8]])
b = np.array([4, 2, 4, 2])
print('a:\n')
print(a)
print('\nb: ', b)
x = spy.linalg.solve_triangular(a, b, lower=True)
print('='*20)
print('x: ', x)
print('\ndot(a, x): ', np.dot(a, x))
.. parsed-literal::
a:
array([[3.0, 2.0, 1.0, 0.0],
[2.0, 1.0, 0.0, 1.0],
[1.0, 0.0, 1.0, 4.0],
[1.0, 2.0, 1.0, 8.0]], dtype=float64)
b: array([4.0, 2.0, 4.0, 2.0], dtype=float64)
====================
x: array([1.333333333333333, -0.6666666666666665, 2.666666666666667, -0.08333333333333337], dtype=float64)
dot(a, x): array([5.333333333333334, 1.916666666666666, 3.666666666666667, 2.0], dtype=float64)
|
