You need a pencil. Pencils are important. And yeah, a paper sheet. A lot of them, actually. And don't forget the course on linear algebra Essence of linear algebra delivered by 3Blue1Brown.
Another useful resources to start are the module on linear algebra and the modules about vectors and matrices in the "Precalculus" section of the Khan Academy, and the legendary Introduction to Linear Algebra (now at its 6th edition) by Gilbert Strang.
This page describes how the instructions vfmadd132ps, vfmadd213ps and vfmadd231ps work.
Wait, what's a floating point value, exactly?
A floating point value is a real number using an integer with a fixed precision, called the significand, scaled by an integer exponent of a fixed base. In computing, floating-point arithmetic (FP) is arithmetic that represents subsets of floating point numbers. From Wikipedia.
Notice that vfmadd132ps, vfmadd213ps and vfmadd231ps are instructions, and not system calls.
An instruction set architecture defines the machine code understood by some processor. An instruction changes the (observable) state of the computer (e.g. changes content of processor registers -including the program counter and the call stack pointer, memory locations in virtual address space, etc...).
A system call, on the other hand, is done by some application program to request services from the operating system kernel. It may correspond to an elementary machine instruction (e.g. SYSENTER or SYSCALL), but the kernel will run a big lot of code before returning to the application program. From Wikipedia.
We can find these instructions in Rust, C and C++.
A Guide to Porting C/C++ to Rust
Before starting, I learnt to use Rust with the Rustlings, a series of exercises which is very helpful to understand Rust and its compiler.
To prepare for this exercise, I implemented the vector and matrix structures and some help methods as print, rows (for vector), shape and is_regular (for matrix).
The goal of the exercise 00 is to implement, for each struct, the methods for addition, subtraction and scaling.
Using generic types implies implementing functions for a specific data type if the operators used are not implemented for that particular type. Traits are used in Rust to implement operator overloading.
For my main function I used a list of colour escape sequences (it's important to note that Rust only considers hexadecimal escape sequences).
If we consider a vector as a movement, adding two vectors is like starting another movement after the first one is ended. The sum of two vectors will then be a vector starting from the origin and reaching the end of the second vector, whose starting point is the tail of the first vector.
In vector addition, each component of the second vector is added to the corresponding component of the first vector.
In the same way, multiplying a vector by a positive number means stretching or squishing it, while multiplying it by a negative number means stretching or squishing it but in the opposite direction.
Multiplying a vector by a scalar means multiplying each of its components by a scalar.
For this exercise 0I had to implement a function calculating a linear combination of an array of vectors and an array of scalars.
To implement the linear_combination function, I used the fused multiply-add function implemented by the mul_add function for the f32 primitive.
Linear combination means adding to a one combination of a scalar and a vector one or more combinations of the same kind. However, if the vector are linearly independent, the set of their linear combinations may span the entire space, while if they are linearly dependent, their span may be limited to a plane (in 3D space), a line (in 2D space, if they line up) or even a point.
However, we can also consider that any of them is outside the span of the others.
From here we can address the technical definiton of a basis:
- the basis of a vector space is a set of linearly independent vectors that span the full space.
Every other vector may be considered as a linear combinattion of the basic vectors.
Linear interpolation is a method of curve fitting (i.e. the process of constructing a curve, or mathematical function, that has the best fit to a series of data points) using linear polynomials to construct new data points within the range of a discrete set of known data points.
A polynomial is a mathematical expression consisting of indeterminantes and coefficients involving only the operations of addition, subtraction, multiplication, and positive-integer powers of variables of degree one, and a linear polynomial is a polynomial of degree one.
For implementing the lerp function for the types f32, Vector<f32> and Matrix<f32> I implemented the Lerp trait and the function lerp derived from it:
pub trait Lerp<V> {
fn lerp(&self, u: V, v: V, t: f32) -> V;
}
impl Lerp<f32> for f32 {
fn lerp(&self, u: f32, v: f32, t: f32) -> f32 {
if t < 0.0 || t > 1.0 {
panic!("t must be comprised between 0 and 1");
}
(1. - t).mul_add(u, t * v)
}
}
impl Lerp<Vector<f32>> for Vector<f32> {
fn lerp(&self, u: Vector<f32>, v: Vector<f32>, t: f32) -> Vector<f32> {
if t < 0.0 || t > 1.0 {
panic!("t must be comprised between 0 and 1");
}
u.mul_add(1.0 - t, &(t * v))
}
}
pub fn lerp<V: Clone>(u: V, v: V, t: f32) -> V
where V: Lerp<V>,
{
u.lerp(u.clone(), v, t)
}I also implemented the mul_add function for Vector<f32> and the std::ops::Mul<Vector<f32>> trait for the f32 primitive to enable multiplication between Vector<f32> and f32:
impl std::ops::Mul<Vector<f32>> for f32 {
type Output = Vector<f32>;
fn mul(self, _rhs: Vector<f32>) -> Vector<f32> {
let mut a: Vector<f32> = Vector::new();
for el in _rhs.values.iter() {
a.values.push(*el * self);
}
a
}
}The same had to be done for Matrix<f32>.
Trait bounds directly inserted in between angle brackets after the impl keyword specify constraints on the type parameters for the entire implementation block. These bounds apply globally to all functions and associated items within the implementation. For example:
impl<K: Clone + std::fmt::Display + std::ops::AddAssign> Vector<K> {
// All methods in this implementation can assume that K implements Clone, Display, and AddAssign.
}In this context, the bounds K: Clone + std::fmt::Display + std::ops::AddAssign mean that any type K used in this implementation of Vector<K> must implement the Clone, Display, and AddAssign traits.
On the other hand, a where clause can be used to specify trait bounds in a more flexible and expressive manner. It allows you to add constraints that are more complex or involve multiple types; it can be applied to individual functions or methods, or to the entire implementation block:
impl<K> Vector<K>
where K: Clone + std::fmt::Display + std::ops::AddAssign,
{
// All methods in this implementation can assume that K implements Clone, Display, and AddAssign.
}In this case, the bounds where K: Clone + std::fmt::Display + std::ops::AddAssign serve the same purpose as placing the bounds directly between the angle brackets of impl<>, but where clauses are preferred when the constraints are more complex or to make the code more readable.
Here are explanations for the dot product taken from Wikipedia and an article from the Khan Academy:
The dot product measures how much two vectors point in the same direction.
From Wiki:
Taxicab norm is the sum of the size of the two catheti of the triangle whose hypotenuse is represented by the vector.
The Euclidean norm is the root square of the sum of the squares of the two catheti of the triangle whose hypotenuse is represented by the vector.
The supremum norm of a vector is the maximum absolute value of its components, or the maximum distance from the origin among those of the components of the vector. It is used in various mathematical and engineering contexts to measure the size of a vector by focusing on its largest component and is particularly useful in optimization and numerical analysis because it provides a simple and intuitive way to assess the magnitude of vectors in terms of their most significant element.
To calculate the cosine of an angle between two vectors cos(θ), I used the formula:
impl AngleCos<f32> for Vector<f32> {
fn angle_cos(&self, u: Vector<f32>, v: Vector<f32>) -> f32 {
u.dot(v.clone()) / (u.norm() * v.norm())
}
}Cosine is the ratio between the size of the adjacent side of the triangle formed by the vector and its hypotenuse (cf. the acronym SOH-CAH-TOA)
The cross product is a binary operation on two vectors in a three-dimensional oriented Euclidean vector space, and is denoted by the symbol ×. Given two linearly independent vectors a and b, the cross product, a × b (read "a cross b"), is a vector that is perpendicular to both a and b, and thus normal to the plane containing them.
The cross product is defined by the formula
a × b = ‖ a ‖ ‖ b ‖ sin( θ ) n
where
θ is the angle between a and b in the plane containing them (hence, it is between 0° and 180°),
‖a‖ and ‖b‖ are the magnitudes of vectors a and b,
n is a unit vector perpendicular to the plane containing a and b, with direction such that the ordered set (a, b, n) is positively oriented.
If the dot product measures how much two vectors point in the same direction, the cross product measures how much two vectors point in different directions.
I implemented the following formula:
impl CrossProduct<f32> for Vector<f32> {
fn cross_product(&self, u: Vector<f32>, v: Vector<f32>) -> Vector<f32> {
if u.rows() != 3 || v.rows() != 3 {
panic!("Vectors must have three dimensions.");
};
let mut cross = Vector::new();
cross.push(u[1].mul_add(v[2], -(u[2] * v[1])));
cross.push(u[2].mul_add(v[0], -(u[0] * v[2])));
cross.push(u[0].mul_add(v[1], -(u[1] * v[0])));
cross
}
}To implement the cross_product function I had to implement the Index trait for Vector<f32>:
impl std::ops::Index<usize> for Vector<f32> {
type Output = f32;
fn index(&self, index: usize) -> &f32 {
&self.values[index]
}
}To talk about matrix multiplication, we need to address the subject of linear transformation, or linear map. A transformation is a function, taking some vector as input and spitting out some vector as output. It is linear if it has two properties:
- all lines must remain lines, without getting curved;
- the origin must remain fixed in place.
In this way, linear transformation is a way to move around space such that gridlines remain parallel and evely spaced, and such that the origin remains fixed.
Following from this, a matrix may be considered as a series of column vectors representing the positions in which the basis (or unit) vectors of the particular dimensional space the matrix belongs to land.
Thus, if we want to know where a particular matrix coordinates take a vector, we multiply each entry of the vector with its corresponding matrix column and add together what we get:
In this way we add the scaled version of our new basis vectors.
In the same way, multiplying two matrices means applying two linear transformations one after the other to obtain a composed linear transformation.
We can see matrix multiplication as follows:
Notice the order of the matrices: the matrix representing the first transformation is put on the right and the one representing the new transformation is put on its left.
Matrix multiplication is associative, meaning that as long as the transformations are applied in the same order, it does not matter whether two matrices are multiplied by themselves first before their product is multiplied by another matrix.
The trace of a square matrix in linear algebra is defined to be the sum of elements on the main diagonal (from the upper left to the lower right) of the matrix. It is only defined for a square matrix.
The trace of an n × n square matrix A is defined as:
where aii denotes the entry on the ith row and ith column of A. The entries of A can be real numbers, complex numbers, or more generally elements of a field F.
The trace is useful for various reasons, among which:
- it is the sum of the eigenvalues of the matrix;
- it has a lot of nice properties such as linearity, invariance by transposition and basis change, and invariance by cyclic permutations; i.e. trace(ABC)=trace(CAB)=trace(BCA) for square matrices A,B,C.
The transpose of a matrix is an operator which flips a matrix over its diagonal; that is, it switches the row and column indices of the matrix A by producing another matrix, often denoted by AT (among other notations).
A matrix is in row echelon form if:
- all rows having only zero entries are at the bottom;
- the leading entry (that is, the left-most nonzero entry) of every nonzero row, called the pivot, is on the right of the leading entry of every row above.
Some texts add the condition that the leading coefficient must be 1 while others require this only in reduced row echelon form.
These two conditions imply that all entries in a column below a leading coefficient are zeros.
The following is an example of a 4 × 5 matrix in row echelon form, but not in reduced row echelon form (see below):
I will make use of Gaussian elimination algorithm to reduce matrices to their row echelon form.
In mathematics, Gaussian elimination, also known as row reduction, is an algorithm for solving systems of linear equations consisting of a sequence of elementary row operations modifying the matrix until the lower left-hand corner of the matrix is filled with zeros. There are three types of elementary row operations:
- swapping two rows,
- multiplying a row by a nonzero number,
- adding a multiple of one row to another row.
Using these operations, a matrix can always be transformed into an upper triangular matrix, and in fact one that is in row echelon form.
Once all of the leading coefficients (the leftmost nonzero entry in each row) are 1, and every column containing a leading coefficient has zeros elsewhere, the matrix is said to be in reduced row echelon form. This final form is unique; in other words, it is independent of the sequence of row operations used. For example, in the following sequence of row operations (where two elementary operations on different rows are done at the first and third steps), the third and fourth matrices are the ones in row echelon form, and the final matrix is the unique reduced row echelon form.
To implement REF and RREF algorithms, I found this answer inspiring.
Let's start from a view on 2-dimensional space.
Since for a transformation to be linear, gridlines must remain parallel and evenly spaced after the transformation, what happens to the single square of area 1x1 formed by the two basis vectors during the transformation applies to every square in the grid.
The scaling factor by which a linear transformation changes any area is called the determinant of that transformation.
For example, a determinant of 3 implies that the associated linear transformation increases the area by a factor of 3.
It may be 0 if the linear transformation squishes everything into a smaller dimension (in a 2-dimensional space onto a line or a point, for example).
A negative determinant implies an inversion of space orientation.
We can compute the determinant of 2x2 matrices as follows:
In 3-dimensional space, the basis vectors form a cube with a volume of 1. Since the determinant gives the factor by which every volume is scaled, we can say that it represents the volume of the parallelepiped in which that cube is turned into.
A volume of 0 means that all space is turned into something with 0 volume, either a flat plane, a line or a single point.
An inversion of space orientation here may be detected with the right hand rule:
If we need our left hand to represent the basis vectors, it means that orientation has changed.
To compute the determinant of 3x3 matrices and higher:
Since matrix multiplication is associative, the following rule applies:
The method I use to compute the determinant is called Laplace expansion.
However, this method has a very high time complexity (more than O(n3)). An implementation which has an inferiour time complexity implies computing the row echelon form of the matrix to obtain an upper triangle form and then calculating the product of the entries of the leading diagonal. Finally, we have to take into account the effects of the elementary row operations on the determinant of a matrix: indeed, after a row operation the determinat of the original matrix may differ from that of the modified matrix. Therefore, we will have to multiply the product of the leading diagonal entries by the multiplier storing the effects implied by the row operations.
fn row_swap(&mut self, a: usize, b: usize) {
let tmp = Vector::capture_row(self.clone(), a);
self[a] = self[b].clone();
self[b] = tmp.get_values();
self.det_multiplier *= -1.;
}
fn add_row_multiple(&mut self, a: usize, b: usize, scalar: f32) {
self[a] = (scalar * Vector::from(self[b].clone()) + Vector::from(self[a].clone())).get_values();
}
<SNIP>
pub fn determinant_by_row_reduction(&self) -> f32 {
let rref = self.row_echelon_form();
let mut product = 1.;
for row in 0..rref.rows {
product *= rref[row][row]
}
return product * rref.det_multiplier;
}From the subject:
If you remember the aside on the exterior algebra from before, you might want to know that n-vectors inside an n-dimensional exterior algebra tend to behave a lot like scalars, and are referred to as "pseudoscalars" for this reason.
The determinant is actually the magnitude of the pseudoscalar (the measure of the n-parallelepiped) which is created by the successive wedge product of all the columns vectors of a square matrix. Read that again, slowly.
The inverse or reciprocal of a matrix is a matrix which on multiplication with the given matrix gives the multiplicative identity.
For a square matrix A, its inverse is A-1, and A · A-1 = A-1· A = I, where I is the identity matrix.
The matrix whose determinant is non-zero and for which the inverse matrix can be calculated is called an invertible matrix.
The inverse of a matrix of size 2x2 may be calculated as follows:
Hint: d is the inverse of a and viceversa.
To calculate the inverse of a matrix of greater size, we need to find the adjoint matrix of our matrix, that is calculated as follows:
Firstly, we compute its matrix of minors:
secondly, we compute its matrix of cofactors:
by following the following scheme:
and finally, we calculate the adjugate (or adjoint) by transposing the elements of the matrix of cofactors:
The rank of a matrix can be defined as the largest order of any non-zero minor in the matrix.
It's also the number of linearly independent vectors in the matrix. That's why a non-zero determinant in a square matrix tells us that all rows (or columns) are linearly independent vectors, so it is "full rank" and its rank equals the number of rows.
Indeed, to calculate it, we can compute the row echelon form of the matrix and count the number of entries in the diagonal (a zero line meaning a vector which is linearly dependent from another).
In 3D computer graphics, a viewing frustum or view frustum is the region of space in the modeled world that may appear on the screen; it is the field of view of a perspective virtual camera system.
The view frustum is typically obtained by taking a geometrical frustum — that is a truncation with parallel planes — of the pyramid of vision, which is the adaptation of (idealized) cone of vision that a camera or eye would have to the rectangular viewports typically used in computer graphics. Some authors use pyramid of vision as a synonym for view frustum itself, i.e. consider it truncated.
The exact shape of this region varies depending on what kind of camera lens is being simulated, but typically it is a frustum of a rectangular pyramid (hence the name). The planes that cut the frustum perpendicular to the viewing direction are called the near plane and the far plane. Objects closer to the camera than the near plane or beyond the far plane are not drawn. Sometimes, the far plane is placed infinitely far away from the camera so all objects within the frustum are drawn regardless of their distance from the camera.
Viewing-frustum culling is the process of removing from the rendering process those objects that lie completely outside the viewing frustum.[6] Rendering these objects would be a waste of resources since they are not directly visible. To make culling fast, it is usually done using bounding volumes surrounding the objects rather than the objects themselves.
We need to calculate a projection matrix for the camera frustum.
Projection matrices are specialized 4x4 matrices designed to transform a 3D point in camera space into its projected counterpart on the canvas. Essentially, when you multiply a 3D point by a projection matrix, you determine its 2D coordinates on the canvas within NDC (Normalized Device Coordinates) space.
For my algorithm I followed this very good explanation by Mads Elvheim on StackOverflow.
The clip matrix is the projection matrix.
The projection function will create a projection matrix from 4 parameters:
- the field of view, the angle of the cone of vision,
- the window size ratio
$\frac{w}{h}$ where w and h are the width and height of the rendering window, - the distance of the near plane,
- the distance of the far plane.
The field of view (FOV) is the angular extent of the observable world that is seen at any given moment. It is typically specified in degrees (that's why I need to convert it to radians in my function).
pub fn projection(fov: f32, ratio: f32, near: f32, far: f32) -> Matrix<f32> {
let f: f32 = 1.0 / f32::tan(fov * (PI / 180.0) * 0.5);
return Matrix::from([
[f * ratio, 0., 0., 0.],
[0., f, 0., 0.],
[0., 0., -(far + near) / (far - near), -1.],
[0., 0., -(2.0 * near * far) / (far - near), 0.]
])
}To test my function, I implemented a function to write an array of floats into a file. To do that, I open a file, loop through the float values and convert each one of them into a string before writing it into the file.
pub fn write_matrix_to_file(&self, filename: &str) {
let mut file = File::create(filename).unwrap();
for v in self.values.iter() {
for (j, n) in v.iter().enumerate() {
let round = (n * 100.).round() / 100.;
let mut tmp = String::from(round.to_string());
if j < v.len() - 1 {
tmp += ", ";
} else {
tmp += "\n";
}
file.write(tmp.to_string().as_bytes()).unwrap();
}
}
}
fn main() {
let P: Matrix<f32> = projection(45.0, 1.2, 0.1, 100.);
P.write_matrix_to_file("../varia/matrix_display/proj");
P.print();
}Complex numbers are numbers made up of a real number and an imaginary part i. i is the imaginary unit.
Complex numbers are used to take the square root of negative numbers.
In Rust, the Complex type has many useful methods we can use.
The dot product between complex vectors is given by the dot product of the first term with the conjugate of the second term, i.e. the number with an equal real part and an imaginary part equal in magnitude but opposite in sign.
To calculate the Euclidean norm or L2 norm of a complex number, we need to find the absolute value, or norm, of each entry of the complex vector:
and the same is also true for L1 norm and infinity norm
Here is some help to work with Complex number type in Rust, and here may be found hints to implement one's own complex number class.
To create tests, I took inspiration from the work of Glagan on GitHub.
Markdown supports mathematical notation.
Passing arguments by reference may be used to avoid copying them in the function scope. At the same time, traits bounds which ask for recursivity won't need a copy of the arguments if they are passed by reference to the argument, as can be seen by the implementation of the cross_product method:
pub trait CrossProduct<K> {
fn cross_product(&self, u: &Vector<K>, v: &Vector<K>) -> Vector<K>;
}
<SNIP>
pub fn cross_product<K: Clone>(u: &Vector<K>, v: &Vector<K>) -> Vector<K>
where Vector<K>: CrossProduct<K> {
u.cross_product(u, v)
}As concerning time and space complexity, the notation used is called big O notation, useful in computer science when analyzing algorithms for efficiency.
We have to count the number of operations, or steps, in relation to a constant n given in this case by the number of dimensions of our vectors and matrices and, in the same way, we have to calculate the space in memory used by our functions in relation to the constant n.
To test, you can run cargo test in the root folder or cargo run in any exercise_* directory.
Results may be compared with those of this complex matrix calculator accepting both real and complex values as entries.
I used assert_eq!() expression as suggested in this interesting discussion on StackOverflow.
I could adapt the approx_eq macro to my containers, but I didn't consider it as necessary at this stage.