MATLAB Introductory Guide
May 14, 2025
Table of Contents
Introduction
disp('Hello, world!')
MATLAB, a portmanteau of Matrix Laboratory, is a special-purpose programming language created by Cleve Moler and developed by MathWorks. Its origins lie in LINPACK, a numerical linear algebra library. Since then, it has evolved into a powerful tool for scientific computation, particularly with the Simulink add-on which is used for modeling dynamic systems.
Scientific computing is utilizing the capabilities of a computer to solve difficult problems. It is an example of numerical methods, where you generally approximate mathematical processes. It is at the intersection of experimentation and theory, and applies principles from both applied mathematics and computer science (it should not be confused with computer science or computer engineering). It is so prominent as a discipline that some scientists choose to spend their careers working in computational science. There are entire institutes and research groups dedicated to this.
While similar tools Mathematica and SymPy are better for symbolic manipulation, MATLAB excels at numerical computation. It is thus used in industry, academia, and national labs across a variety of disciplines. Electrical, aerospace, and mechanical engineers are among its most common users.
Some applications of MATLAB include:
- Structural dynamics
- Controls engineering
- Signal proccessing
- AI/ML
- Biology, such as in gel electrophoresis and neuroimaging
- Climate modeling
- Exoplanet imaging
- Temperature Programmed Reaction Spectroscopy
- Most of your math courses can be done in MATLAB
MATLAB, however, is nowhere near the fastest language out there. Languages like R, Python, and Fortran can be just as fast, if not faster, and sometimes more common in certain industries. C/C++ can be orders of magnitude faster. That being said, MATLAB remains as one of those most prominent environments for prototyping—testing out and visualizing your early ideas. At the end of the day, the best tool is subjective, and dependent on your field. Having an arsenal of programming languages will only give you more opportunities, and MATLAB is a great one to start with for learning scientific computing.
This guide was created in MATLAB version R2024B using an interactive Live Script document, which is a great feature for documenting your learning.
1. Using Variables and Functions
Open up MATLAB and navigate to the Command Window. This is kind of like the terminal of MATLAB. MATLAB can perform calculations like arithmetic:
5 + 5
Notice how this answer appears in the Command Window, but also in Workspace and Command History. You can clear Command Window by inputting clc and Workspace by clear:
One of the big topics in programming is being able to generalize functions to make processes easily repeatable and scalable. The first step is to assign variables. In MATLAB, the equals sign (=) is used to assign values to variables. It does not mean the same thing that it does in mathematics as you can see:
format debug
x = 5 + 5
With x now assigned to 10, notice how there is a given structure address for this action when we use format debug. The assignment of a value to a variable allocates a spot in memory. The way data are stored in MATLAB makes it optimal for its original purpose which was for linear algebra. The information given by format debug is not always needed, so for now we can use format short to revert the display to how it is by default.
Let's try calling x:
format short
x
As done before with x, will assign the value 8 to a new variable y. Then, we will add 3 to y and store it as x, overwriting the previous action on line 4:
y = 8
x = y + 3
With two variables that can be called at any time, we can easily perform calculations with whatever the values may be at that point in the program. Let's just try multiplying them together:
m = x * y
MATLAB has a lot of built-in functions that are only available through add-ons or libraries in a language like C++. Some of these can be quite complex, so let's look at some basic ones you are familiar with, the trigonometric functions. Here, we assign the output of the cosine function when a value in degrees is inputted:
f = cosd(180)
We can also plug variables into these functions and manipulate the variables within them. We will call the previous function explicitly into a different trigonometric function, inverse tangent in degrees:
t = atand(cosd(180))
Let's look at a very basic few lines of code that finds and stores the prime factors of 127, divides it by 6, and then rounds the value:
p = factor(127)
d = p / 6
r = round(d)
It works, but it's not self-documenting, and you may run out of letters to use for variables later on. Self-documenting code is code that tells the reader a little bit about how the program works by simply existing. Let's rewrite this code with self-documenting variable names:
factored_value = factor(127)
division_step = factored_value / 6
final_rounding = round(division_step)
It might not look like much now, but in a file with potentially thousands of lines of code, a little bit of readability goes a long way, especially when working with peers. Some general rules for naming variables is to make sure they begin with a letter, cannot contain punctuation, cannot contain spaces, cannot be named as MATLAB Keywords, and should not be named as MATLAB built-in functions. Doing so will overwrite the function as a variable, rendering the function unuseable. Investigate this on your own to see what happens.
You can test to see if something is available as a variable by using exist:
exist global
It looks like global is returning 5, so it is reserved for a function should we decide to call it. Let's look it up to verify:
help global
Indeed, global is already a built-in function.
Let's try a different variable name:
exist universe
universe returns 0, so it is fair game.
There are different types of variables in MATLAB, which can be seen in the Workspace tab. You have to right click the tab and select Class for it to show an additional column of data. Let's try this out:
x = 4^2
This means that the variable x is assigned the value 16 using the default data class "double." This is known as a double precision floating point number. You can typecast values to different types of data classes:
x_single = single(4^2)
Notice how the new class here is single.
By right clicking on a variable in Workspace, you can open up the Variable Browser to view and edit your variables:
Let's manually add in new values and see what happens:
Now, let's call x:
This also results in a new address for x.
This can also be done in the Command Window. This is known as a row vector:
basic_row_vector = [16, 2, 5, 47]
Strings are a data class that represent a sequence of characters. They can be stored in row vectors too:
stellar_remnants = ["neutron_star", "black_hole", "white_dwarf"]
There will be more on arrays very soon. Before that, let's try writing a very basic program as script in MATLAB that calculates the Arrhenius equation. This is one of the most important equations in chemistry, which relates temperature to the reaction rate:
Notice how this is a script file, which can be used to execute more than one line of code multiple times when saved.
Try doing this for an equation in a subject you are studying. Using MATLAB's documentation is essential to learning the language. I highly recommend using it here. Try picking an equation that will require you to use a new built-in function.
2. What Is an Array?
Every computer and program uses arrays to some capacity. An array is a data structure that consists of a collection of elements. The mathematic term for an array is a matrix. Algebra is a branch of mathematics that is all about analyzing and then manipulating sets of things, which are formally known as algebraic structures. Thus, linear algebra is about the analysis of linear equations and linear systems, which manifest in vector spaces and matrices. Numerical linear algebra is one of the most prominent topics in mathematics today, with significant applications in data science, electrical engineering, quantum mechanics, and more. Perhaps most recognizable is the fact that AI is entirely dependent on linear algebra. If you haven't already taken a linear algebra course, chances are, you will in the next few years.
A scalar is just a 1 x 1 matrix.
Row vectors look like this:
And column vectors like this:
These are among the first and most common types of arrays.
This is how you create a row vector:
row_vector_sample = [1, 0, 1, 0]
This is how you create a column vector:
column_vector_sample = [1; 0; 1; 0]
Let's look at this data set, which represents data taken for Ohm's law, which relates electric current to voltage:
Current (Amps) | Voltage (V) |
0.00 | 0.00 |
1.00 | 2.31 |
2.00 | 3.94 |
3.00 | 6.23 |
4.00 | 7.28 |
5.00 | 9.67 |
First, let's assign these values to two vectors:
current = [0.00, 1.00, 2.00, 3.00, 4.00, 5.00]
voltage = [0.00, 2.31, 3.94, 6.23, 7.28, 9.67]
.Now, we can plot them:
plot(current, voltage, '-or')
xlabel('Current (Amps)')
ylabel('Voltage (V)')
title('Voltage vs. Current of a Resistor')
Technically, you can use either row or column vectors for the same things. It just comes down to standards for data representation and what works for you.
Now, let's introduce the colon operator (:). This is one of the most useful operators available. It lets you create vectors for a certain number of iterations.
We have the following row vector:
no_colon_row_vector = [0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20]
It works fine, but what happens when we have a data set much larger than this one? Inputting that manually was also pretty tedious.
Take a look at the first usage of the colon operator:
x = j:i:k creates a vector x with increments of i that spans from j to k
Now, let's apply it:
colon_row_vector = [0:2:20]
% The vector starts at 0, ends at 20, and has 2 spacing between each element.
We are able to create the same row vector much faster. With this much power, you can create very large vectors. In cases like the one below, it's best to suppress the output with a semicolon (;):
giant_vector = [1:5:100];
The vector is still stored, just not displayed.
There are limits to colon notation. Let's say we have i set to a value that will occur before the upper bound is reached:
1800:50:2025
Notice how the last value displayed is 2000, rather than the 2025 we set. There is no way to get to 2025 from 1800 if we are adding by increments of 50. Colon notation is for vectors with equal increments, also known as evenly spaced vectors.
What happens if we set i to a value that is in the opposite direction of the vector?
2019:-10:2049
It is impossible to get from a lower bound of 2019 to an upper bound of 2049 by going down 10 increments each time. Thus, no information is stored. This returns an empty vector.
Try to think of some interesting things you can do with row vectors. Here's a basic program that adds two vectors data_set1 and data_set2, and only returns even integers:
% Two row vectors
data_set1 = [2 3 4 5 6]
data_set2 = [1 5 6 11 12]
% Add them together
data_set_combined = data_set1 + data_set2;
% The modulo function calculates the remainder of division.
% Even numbers return a remainder of 0 when divided by 2.
% This returns values from data_set_combined that meet this condition.
evens = data_set_combined(mod(data_set_combined, 2) == 0)
Let's learn about the linspace function:
help linspace
In other words, this function is helpful when you need a certain number of points, whereas colon notation is for step size.
Here's a brief review question. Which of the following will NOT output a vector of 10 elements? Commit to an answer and then test the following code yourself:
optionA = [1:1:10];
optionB = [0, 100, 10];
optionC = linspace(0, pi, 10);
optionD = linspace(7.5, 7.51, 10);
Now, we can learn about matrices. A matrix is a rectangular array of real numbers. Let m = rows and let n = columns for this m by n matrix A:
Matrix M,
,can be created like this:
M = [1, 2, 3 ; 4, 5, 6 ; 7, 8, 9] % In math, matrices are assigned uppercase letters and vectors lowercase letters.
We can also generate it using colon notation, so we don't have to do it entry-by-entry:
M = [1:3; 4:6; 7:9]
3. Plotting
Before we learn more about arrays, let's provide an overview on how to plot.
Here are some plot elements:
Argument | Line |
- | Solid line |
-- | Dashed line |
: | Dotted line |
-. | Dash-dotted line |
Argument | Symbol |
o | Circle |
+ | Plus |
* | Asterisk |
. | Point |
x | Cross |
square | Square |
^ | Upwards triangle |
Argument | Color |
y | yellow |
m | magenta |
c | cyan |
r | red |
g | green |
b | blue |
w | white |
k | black |
Let's try making a sin graph:
% First, we must establish a set of x-values. We are going from -pi to pi.
x = -pi:0.01:pi;
% Here is the function we are graphing.
y1 = sin(x);
% Now, we can plot it. Remember to follow proper graphing standards.
plot(x, y1, 'b-') % b denotes the line being blue and - denotes a solid line
xlabel('x (radians)')
ylabel('sin(x)')
title('Plot of y1 = sin(x)')
We can of course have multiple sin waves on the same axis. Let's define another and then graph it in green:
y2 = sin(x + pi);
plot(x, y1, 'b-', x, y2, 'g-')
xlabel('x (radians)')
ylabel('sin(x)')
title('Plot of Sin Waves')
Let's set a legend:
legend('sin(x)', 'sin(x + pi)')
There are other types of plots, like this 3-D graph:
% Define the range of x and y values
x = -2*pi:0.1:2*pi;
y = -2*pi:0.1:2*pi;
% Create a grid of x-y coordinates
[x, y] = meshgrid(x,y);
% Here is our function
z = sin(x).*cos(y);
% Create the 3D surface plot
surf(x, y, z)
% Label axes
xlabel('x')
ylabel('y')
zlabel('z = sin(x) * cos(y)')
title('3D Graph of sin(x) * cos(y)')
% Add shading and a color bar
shading interp
colorbar
4. Working with Arrays
Let's learn how to traverse a matrix. First, let's populate an arbitrary 5 x 5 matrix using rand:
help rand
A = rand(5,5)
You can easily access the number of rows and columns by using size:
size(A)
We can call 
:
:A(1,2)
You can import a matrix from other programs using the file format Comma-separated values (.csv). Simply use Import Data with Output Type set to Numerical Matrix.
As we index, the colon operator further proves its value. Take a look at the syntax:
Let's test out some of these.
Here's calling column 3:
A(:,3)
Here's calling row 2:
A(2,:)
Here's calling rows 4 to 5 in column 1:
A(4:5,1)
Here's calling a small portion of A:
A(1:2, 1:2)
Let's save column 5 as a brand new variable:
c5 = A(:,5)
Let's save row 4 to a brand new variable:
r4 = A(4,:)
Let's convert row 4 to a column vector:
r4(:)
Let's call some very specific elements: columns 4 and 5 of rows 2 to 3:
A([2:3], [4, 5])
Performing operations on matrices is fundamental to modern day science and mathematics.
We can perform matrix addition for the following matrices A and B:
A = [1:3; 4:6]
B = [6:-1:4; 3:-1:1]
C_add = A + B
We can calculate the Hadamard product. This multiplies element in A by the corresponding element in B:
C_Hadamard = A.*B
We can use the array power operator (.^) which raises each element of A to the power of the corresponding element in B:
C_power = A.^B
And the most important is matrix multiplication. For this, we have to pick the proper sized matrix. The number of rows and columns should be inverted. A is fine, but we will pick a new second matrix:
Now, let's test it:
S = [1 0; 0 1; 1 -1];
AS = A*S
Here are two functions that are used for generating matrices.
The first, zeros, populates a matrix with 0s where the syntax is zeros(m,n):
Z_0s = zeros(2,3)
The second, ones, populates a matrix with 1s where the function is ones(m,n):
Z_1s = ones(2,3)
An identity matrix is a special type of matrix where there is a staircase of ones and zeros everywhere else. Multiplying any matrix by the identity matrix is analogous to multiplying by 1 because you will get the original matrix. When the row index is equal to the column index, the element is equal to 0. Everywhere else, it is equal to 1. An identity matrix can be called like this:
I5 = eye(5)
Let's take a look at a matrix where the row and column indices are different:
I23 = eye(4,8)
Another useful operator is the transpose operator (' or transpose). In the transpose of a column vector, the ith row becomes the ith column. In the transpose of a row vector, the kth column becomes the kth row. In other words, a column vector becomes a row vector and vice-versa.
MATLAB prefers row vectors by default:
x = (-4:2:4)
We can transpose this to turn it into a column vector:
x = (-4:2:4)'
Alternatively,
x = transpose(-4:2:4)
Similarly, transpose interchanges the row and column index for each element in a matrix. Let's generate a matrix using magic to showcase this:
X = magic(4)
X_new = transpose(X)
We can also use gallery to call from a library of test matrices, such as dramadah:
X_dramadah = gallery('dramadah', 4)
X_new = transpose(X_dramadah)
Let's try again with a unique matrix known as a Hankel matrix. This can be called using ipjfact as a parameter:
X_Hankel = gallery('ipjfact', 6)
X_new = transpose(X_Hankel)
Transposing matrices has some niche yet interesting applications.
- The transpose of an orthogonal matrix is equal to its inverse
- Transposition ciphers are a part of cryptography
- Transposing can move a matrix from one side of a dot product to the other
- Optimizing data for easier computation
- Signal processing in MIMO systems, the basis of wireless communication
We can also add or remove elements in a matrix.
For the following,
X = [1, 2; 3, 4]
We can add in a third row:
X(3,:) = [5, 6]
We can also add in a third column:
X(:,3) = [7; 8; 9]
We can swap out a specific element. Let's replace 9 with 10:.
X(3,3) = 10
We can remove row 2:
X(2,:) = []
We can remove column 1:
X(:, 1) = []
We can remove a single element only from a vector, not a matrix. Let's get rid of row 2 so we have a vector, and then remove the scalar 7:
X(2,:) = []
X(2) = []
All that remains is the scalar 2. We could have of course been left with a vector had there been more entries.
Here's another example. Let's start with a matrix of entirely zeros:
Z = zeros(5)
Let's eliminate the bottom 3 rows:
Z(1:3,:) = []
Let's append ones to columns 4 and 5:
Z(:,4:5) = ones(2,2)
Remember to be careful with the dimensions when you are adding or removing things from a matrix. You can only add rows with the same number of columns as a matrix and vice-versa. You can't just delete scalars from the middle of a matrix either.
Let's try a little bit of concatenation. We can combine row or column vectors to form a matrix:
r1 = [1, 2, 3];
r2 = [4, 5, 6];
A = [r1; r2]
We can also form a block matrix. A block matrix is a matrix that is broken into sub-matrices as such:
A = [1:3; 4:6; 7:9]; % We have the following matrix.
% We can then generate different sub-matrices and append them onto a new matrix B.
B = [A, zeros(3); ones(3), eye(3)]
A new function for matrixes is diag. Let's define another identity matrix and call it with diag:
I = eye(5);
diag(I)
Recall how an identity matrix has a "staircase" of ones. As you can see, the function returned only ones.
We can also use diag to create a matrix with a certain diagonal. Let's define one matrix and then create another using the first as its diagonal:
a = 0:2:8;
diag(a)
Let's look at the upper and lower triangle entities functions, triu and tril respectively. For fun, we will demonstrate using more matrices from gallery:
A = gallery('fiedler', 5, 1)
L = triu(A)
L = tril(A)
Now, let's do it with a different matrix:
B = gallery('gcdmat', 5, 1)
L = triu(B)
L = tril(B)
The three functions diag, triu, and tril also have a second parameter, k. It can be used to set information regarding where the location of elements in a matrix come from:
v = ones(1,5); % A row vector of five 1s
I_b2 = diag(v,-2) % Identity matrix where the 1s are two diagonals shifted below main
I_b1 = diag(v,-1) % Identity matrix where the 1s are one diagonal shifted below main
I0 = diag(v) % Identity matrix where the 1s are on the main diagonal
I_a1 = diag(v,1) % Identity matrix where the 1s are one diagonal shifted above main
I_a2 = diag(v,2) % Identity matrix where the 1s are two diagonals shifted above main
What do you think it does for triu and tril?
5. Logical Operations
Python has data classes like str, int, and float. Similarly, there are many different types of data in MATLAB. These data classes are of course stored as matrices and scalars. As per usual, the documentation is by far the best way to learn this content.
You already have likely used these:
Type | Purpose |
double | default numeric |
string, char | sequence of characters |
While this section will focus mostly on one data class, let's briefly take a look at two other fundamental classes.
A cell array can contain a combination of different types of elements:
c = {'A B C', [1, 2, 3], 1970}
timetable is used when there is a time associated with each row in a table:
time = datetime(2025,1,1) + days(0:2);
temp = [59.1; 57.8; 58.4];
weather_balloon = timetable(time', temp, 'VariableNames', {'High Temperature (°C)'})
Before we start with the main focus of this section, now is a good time to summarize the terminology for arrays. Recall that an 
matrix has 
and 
. All matrices are arrays, but not all arrays are matrices.
matrix has and . All matrices are arrays, but not all arrays are matrices.Data type | Size |
Scalar |  |
Row vector |  |
Nonscalar row vector |  |
Column vector |  |
Nonscalar column vector |  |
Matrix |  |
Nonscalar matrix |  |
Nonvector matrix |  |
Now, let's talk about the logical data type. It stores one of two values: either true or false. This is more commonly known as the Boolean, named after English mathematician George Boole. His system is the basis for basically every digital computer. In MATLAB, this data type can be stored as a full or sparse matrix. Logical values take up very little storage: only 8 bits. Thinking about it as a variable, the assigned value is stored as an 8-bit string in memory. It can be used to logical operations, relational operations, test the state of variables, and index arrays. Logical is also a function that converts numeric values to logical ones. Let's look at an example:
zero = 0;
one = 1;
T = logical(one)
F = logical(zero)
In Workspace, notice how both the number of Bytes and Class have changed. These functions were not given the proper parameters for a reason. You are supposed to input true or false into that function. We used one and zero because 1 represents true and 0 represents false in coding. When it comes to boolean operations, also known as boolean operations, there are logical operators we can use to perform them. Here's a basic overview.
One of these (~) is used to find the logical NOT (more on this can be found below):
~T
This of course outputs 0 because if something is not true it must be false.
Now, let's try a relational operation (==):
5 == 4
This is of course false and outputs 0 as a result. If we correct this expression,
5 == 5
, we get an output of 1 because it is now true. We can repeat this by storing it as a self-documenting variable like test. If you want to call this later in a program, you can do so easily:
test = 5 == 5;
test
The test outputs true. This is a logical scalar.
Now let's try testing the state of variables using isa. T is currently stored as the logical data type, so let's see what happens if we test if T is an instance of (being stored as) another structure, the cell:
isa(T, "cell")
T is of course not stored as a cell, so it returns 0. Let's correct it to output 1:
isa(T, "logical")
It now returns true because T indeed is stored as the logical data type.
Let's try indexing an array. We have the following array:
x = 1:1:10
% Now, let's check where the array is equal to 5.
x == 5
It is only equal to 5 at the fifth entry in x. Thus, it only returns 1 at the fifth entry in the array being used to compare. It is true that the fifth element is 5. It is false for every other number, which is why every other element is 0.
Let's look at another example while also practicing some of our array generation skills. Let's create the following array:
x = zeros(1,10); % Create a vector of zeroes
t = 2*ones(1,5); % Create a vector of 2s t
x(1, 2:2:10) = t.^([0:1:4]) % Raise each element in t to values from 0 to 4 with 1 spacing and append those values to even positions
l_vector = logical(x)
l_vector stores a 1 where there are nonzero values from x and a 0 otherwise. This is known as a logical vector.
Logical scalars are great for a single thing, like one condition, whereas logical vectors are better for analyzing a set of data.
Logical matrices are of course a thing too. Let's take a random matrix and have logical tell us where the odd values are:
A = randi(50, 5, 5)
L = logical(mod(A,2))
Now, let's provide a formal overview of the most important logical operators in MATLAB. You may have seen this in an introductory computer science class as truth tables.
An important logical operator in MATLAB is NOT (~). It flips the truth value of a logical variable. Let's try this with a matrix:
A = eye(3)
B = ~A
Notice how the values swapped. This new array B is called the logical negation of A. Let's generalize the logical NOT operator:
Logical value of A | Truth value of A | → | Logical value of ~A | Truth value of ~A |
0 | false | → | 1 | true |
1 | true | → | 0 | false |
The NOT operator is great for filtering data and inverting conditions. It can save you from poor code (i.e. having to use too many if statements).
Let's move on to the AND operator (&). Here is the truth table:
A | B | A AND B |
T | T | T |
T | F | F |
F | T | F |
F | F | F |
Both A and B must be true for the condition to be true.
Let's test this in MATLAB using logical scalars:
A = true
B = false
A&B
B&A
What should the following output?
A = true;
B = true;
A&B;
We can also just create the entire truth table using logical vectors:
% Create column A and column B
A = logical([1; 1; 0; 0])
B = logical([1; 0; 1; 0])
% Test and create A and B
A&B
% Combine the vectors to form the truth table matrix
ANDTruthTable = [A, B, A&B]
This matches the truth table.
Let's look at it with two matrices:
A = zeros(6,6);
oneA = ones(1,6);
A(1:2:5,:) = [2*oneA; 4*oneA; 8*oneA]
LA = logical(A)
B = zeros(6,6);
oneB = ones(6,1);
B(:,2:2:6) = [2*(oneB), 4*(oneB), 8*(oneB)]
LB = logical(B)
Now let's test them:
LA&LB % This tests where both LA and LB are true
The first element of LA is 1. The first element of LB is 0. Therefore, this outputted 0. If you need practice with truth tables, try tracing each of the 36 comparisons.
We can of course combine the previous two operators, NOT and AND:
(~LA)&LB % This tests where only LB is true
LA&(~LB) % This tests where only LA is true
(~LA)&(~LB) % This tests where both LA and LB are false
The AND operator can thus test for multiple conditions simultaneously.
Let's move on to the OR operator (|). Here is the truth table:
A | B | A OR B |
T | T | T |
T | F | T |
F | T | T |
F | F | F |
Is this an inclusive or exclusive operator?
Let's test this in MATLAB using logical scalars:
A = true
B = false
A|B
What should the following output?
A = false;
B = false;
A|B;
We can also just create the entire truth table using logical vectors:
% Create column A and column B
A = logical([1; 1; 0; 0])
B = logical([1; 0; 1; 0])
% Test and create A or B
A|B
% Combine the vectors to form the truth table matrix
ORTruthTable = [A, B, A|B]
This matches the truth table.
Let's look at it with two matrices:
A = zeros(6,6);
oneA = ones(1,6);
A(1:2:5,:) = [2*oneA; 4*oneA; 8*oneA]
LA = logical(A)
B = zeros(6,6);
oneB = ones(6,1);
B(:,2:2:6) = [2*(oneB), 4*(oneB), 8*(oneB)]
LB = logical(B)
Now let's test them:
LA|LB % This tests where LA or LB is true
If there is a 1 in either LA or LB, it outputs 1. It only outputs 0 when there is no 1 in either LA or LB.
Let's combine operators:
LA|(~LB) % This tests where LA is true or LB is false.
Note that all of the operators above can also be called in function form:
Operator | Function |
~A | not(A) |
A&B | and(A,B) |
A|B | or(A,B) |
In terms of number of operands, NOT is unary; the AND and OR operators are binary.
Now, we're going to go little deeper into logic. MATLAB has many operators that assess things on an element-wise basis. That means they have to pair the corresponding elements in two arrays. Let's look at the compatible sizes for binary logical operators.
I. Two inputs have identical sizes:
II. One input is a scalar:
III. One input is an matrix and the other is an 
column vector. In other words, the column vector has the same number of rows as the matrix:
matrix and the other is an column vector. In other words, the column vector has the same number of rows as the matrix:
IV. One input is an matrix and the other is a 
row vector:
matrix and the other is a row vector:
V. One input is an column vector and the other is a row vector:
column vector and the other is a row vector:
This means there are incompatible sizes for binary logical operators. I will include those here for reference.
I. A and B are two nonvector matrices with two different row sizes:
II. A and B are two nonvector matrices with two different column sizes:
III. Two nonscalar row vectors A and B are of :
IV. Two nonscalar column vectors A and B are of different sizes:
Earlier, we discussed the OR operator (inclusive). Now, let's discuss the XOR operator (exclusive). Here is the truth table:
A | B | A XOR B |
T | T | F |
T | F | T |
F | T | T |
F | F | F |
XOR is only met if exactly one is true. It is looking for one of A or B, but not both. XOR is called as a function.
Let's test this in MATLAB using logical scalars:
A = true;
B = false;
xor(A,B)
Review: What should the following output?
A = true;
B = true;
or(A,B);
Now, what will this output?
xor(A,B);
We can also just create the entire truth table using logical vectors:
% Create column A and column B
A = logical([1; 1; 0; 0])
B = logical([1; 0; 1; 0])
% Test and create A or B
xor(A,B)
% Combine the vectors to form the truth table matrix
XORTruthTable = [A, B, xor(A,B)]
This matches the truth table.
Let's look at it with two matrices:
A = zeros(6,6);
oneA = ones(1,6);
A(1:2:5,:) = [2*oneA; 4*oneA; 8*oneA]
LA = logical(A)
B = zeros(6,6);
oneB = ones(6,1);
B(:,2:2:6) = [2*(oneB), 4*(oneB), 8*(oneB)]
LB = logical(B)
Now let's test them:
xor(LA,LB) % This tests that out of the same element number in LA and LB, only one of them has a nonzero value.
XOR is helpful because it is both communative and associative. If you XOR something with a value, you can XOR it again with the same value to get the original back. Let's look at a very basic application of it in cryptography. Let's look at the string "hi" being encrypted and decrypted:
message = [0, 1, 0, 0, 1, 0, 0, 0; 0, 1, 1, 0, 1, 0, 0, 1]
keystream = randi([0, 1], 2, 8)
We now have our message and the keystream, a sequence of random random characters that will be used to produce our encrypted message.
Let's perform the operation:
ciphertext = xor(message, keystream)
Now, let's decrypt:
decrypted = xor(ciphertext, keystream)
Let's subtract decrypted from message. If they are identical, it should return a matrix of entirely zeros:
test = message - decrypted
The cipher worked successfully.
There are also more formal ways of combining logical operators. An example of this can be found in the NAND operator. Here is the truth table that compares AND and NAND:
A | B | A AND B | A NAND B |
T | T | T | F |
T | F | F | T |
F | T | F | T |
F | F | F | T |
NAND is only true when AND is false.
NAND is created in MATLAB by combining other logical operators NOT and AND as such:
~(A&B)
Let's create the NAND truth table using logical vectors:
% Create column A and column B
A = logical([1; 1; 0; 0])
B = logical([1; 0; 1; 0])
% Test and create A or B
~(A&B)
% Combine the vectors to form the truth table matrix
NANDTruthTable = [A, B, ~(A&B)]
We can also combine NOT and OR to create the NOR operator. Here is the truth table that compares OR and NOR:
A | B | A OR B | A NOR B |
T | T | T | F |
T | F | T | F |
F | T | T | F |
F | F | F | T |
Let's create the NOR truth table using logical vectors:
% Create column A and column B
A = logical([1; 1; 0; 0])
B = logical([1; 0; 1; 0])
% Test and create A or B
~(A|B)
% Combine the vectors to form the truth table matrix
NORTruthTable = [A, B, ~(A|B)]
Lastly, we can combine NOT and XOR to create the XNOR operator. Here is the truth table that compares XOR and NXOR:
A | B | A XOR B | A NXOR B |
T | T | F | T |
T | F | T | F |
F | T | T | F |
F | F | F | T |
Let's create the XNOR truth table using logical vectors:
% Create column A and column B
A = logical([1; 1; 0; 0])
B = logical([1; 0; 1; 0])
% Test and create A or B
~xor(A,B)
% Combine the vectors to form the truth table matrix
NORTruthTable = [A, B, ~xor(A,B)]
Before we continue, here's a summary:
Operator | Description |
NOT | Returns true if the input is false |
AND | Returns true if all inputs are true |
OR | Returns true if at least one input is true |
XOR | Returns true if an odd number of inputs are true |
NAND | Returns true if not all inputs are true |
NOR | Returns true only if all inputs are false |
XNOR | Returns true if an even number of inputs are true |
These are the standard symbols for logic gates. They are used in schematics for electronics:
And combined, here are all of the logical operators in a truth table:
NOT | NOT | AND | OR | XOR | NAND | NOR | XNOR | ||
A | B | ~A | ~B | A&B | A|B | xor(A,B) | ~(A&B) | ~(A|B) | ~xor(A,B) |
T | T | F | F | T | T | F | F | F | T |
T | F | F | T | F | T | T | T | F | F |
F | T | T | F | F | T | T | T | F | F |
F | F | T | T | F | F | F | T | T | T |
And here they are represented as a logical matrix:
[A, B, ~A, ~B, A&B, A|B, xor(A,B), ~(A&B), ~(A|B), ~xor(A,B)]
Another English mathematician, Augustus De Morgan, created De Morgan's laws. These are foundational to Boolean algebra and logic as a whole:
not(A or B) = (not A) and (not B)
not(A and B) = (not A) or (not B)
These identities are helpful for simplifying expressions. They are an example of equivalent propositions.
Together, these logical operators combine into a framework that is important for computer science on both a theoretical and applied level. Logic gates, for instance, are fundamental to computer architecture.
6. Relational Operators and Determining State
Relational operations compare elements in two arrays to see how the values relate to one another. They will output logical data. To keep things simple, we will first start with numerical values. Here is an overview of the relational operators:
Symbol | Description | Operand type | Function Form |
== | Equal to | Binary | eq |
~= | Not equal to | Binary | ne |
< | Less than | Binary | lt |
<= | Less than or equal to | Binary | le |
> | Greater than | Binary | gt |
>= | Greater than or equal to | Binary | ge |
Let's try these out. We will use a new function here:
fun = @(x) x.^2; % Create the function
% Let's evaluate the integral of it from 0 to 1 and see if it is equal to 5
integral(fun,0,1) == 5
It is not equal to 5.
Let's try it again with the actual result:
integral(fun,0,1) == 1/3
Now that it is correct, it returns 1 as expected.
Here is what it looks like when you test equality for row vectors:
a = 0:10
b = linspace(0,10,11)
a == b
Indeed, these two row vectors are equivalent.
Here we have two matrices:
A = [2, 4; 6, 8]
B = [2, 1; 3, 2^3]
A == B % This shows us where A and B are equal
A ~= B % This shows us where A and B are not equal
We can combine relational and logical operations:
A == 3 | B == 3 % This shows us where either A or B are equal to 3
We will now move to the inequality operators:
A = [1, 3; 5, 7]
B = [1*5, 3*5; 5*5, 7*5]
A > B % This shows us where A is greater than B
A >= 3 & B < 25 % This shows us where A is greater than or equal to 3 and B is less than 25
The isa function also allows us to do this for characters. For tf = isa(A,dataclass), it returns true if A is an instance of dataclass or a subclass of dataclass:
Let's try this out:
randomstring = "This is a string"
tf1 = isa(randomstring, 'char')
tf2 = isa(randomstring, 'string')
tf1 returns 0 because randomstring is not an instance of char. tf2 returns 1 because randomstring is an instance of string.
Let's use isa to test when something is a subclass of a data type:
randomdouble = 5;
whos randomdouble
tf = isa(randomdouble, 'numeric')
double is a subclass of numeric, so isa returns true.
The isa function is part of a larger set of functions called is* functions. Let's look at an example of another:
This tests the version of MATLAB you have. I am not on macOS, so it returns 0. I am however, on PC, so it returns 1 there. This could be helpful for writing different versions of programs depending on the platform someone is using.
There are many more of these functions and I recommend checking out the documentation to learn more.
7. Logical Indexing
Logical indexing is when you select certain elements from an array that satisfy certain conditions.
Let's work with a random matrix A. We will find the large elements. "Large" will be arbitrarily defined as anything larger than 15:
A = randi(20, 6, 6)
largeValues = A > 15 % This assigns 1 to the large elements in the matrix
A(largeValues) % This lists the large values that were found.
We can now delete those values if desired:
A(largeValues) = 0
With this, let's try to filter even more. We can use the any function to determine if any array elements are nonzero.
Let's perform another test for values that are small. "Small" will be arbitrarily defined as anything less than or equal to 9. We will also ignore the elements that have already been assigned 0:
smallValues = A <= 9 & A ~= 0
test2 = any(smallValues)
This gives us the columns that have at least one nonzero value in them:
On the other hand, the all function returns true if every element is nonzero:
B = [0 0 3;0 3 3;0 0 3]
all(B)
Only the third column returns a 1 because it is the only one comprised of entirely nonzero elements.
You may have noticed that the any and all functions output an array that is a different size from the input array:
Input size | Output size |
|  |
| |
| |
You might come across special values Inf and NaN while working with matrices. We can produce them manually as such:
1/0
-1/0
syms x
double(limit(1/x, x, 0, 'right'))
They can also be operated with:
Inf + 5
Here are some ways to produce NaN (Not a Number):
0/0
Inf/Inf
NaN appears as a placeholder for missing or invalid data that you may have imported from something like a spreadsheet or lab instrument. NaN data appear as holes in graphs:
x = linspace(0, 10, 100);
y = sin(x);
% Introduce a hole
y(40:50) = NaN;
plot(x, y, 'b-', 'LineWidth', 2)
xlabel('x')
ylabel('y')
title('sin(x) Function')
grid on
Let's say we have a matrix containing some Inf and NaN values:
A = [1 2 Inf; 4 NaN 6]
search1 = A > 10^9
search2 = A ~= A
Inf was returned as true when calling for values greater than 10^9.
NaN was returned as true when calling for values in A that are not equal to themselves. This is a little bit of a weird concept:
NaN ≠ NaN
Thankfully, we have the more streamlined isinf and isnan functions to find problematic elements in a matrix.
8. Control Flow
There are going to be times where you want the individual instructions of your program to be executed in a certain order. This is referred to as control flow. Some of the things you might want to do include the repetition of a certain statement until some condition is met, continuing the program at a different statement, or calling certain functions only if a specific condition is true. Conditionals are the constructs that guide your program down certain paths. The four used in MATLAB are for, if, while, and switch.
Let's start with the for loop. This is the general form in pseudocode:
for index = array
statements
end
It will loop over the condition, which is often a numerical array, executing the statements once for each element in it. Here is an example:
for i = [1 2 3 4]
fprintf('The value of index i is %d \n', i)
end
Let's try another example where we start adding to count at 5:
count = 5;
for x = 1:8
count = count + 1;
fprintf('The count is %d \n', count)
end
Here, we approximate the area under a curve using a left Riemann sum. Here is a visualization:
Here's the code:
a = 0; % Start of interval
b = 2; % End of interval
n = 10; % Number of subintervals (rectangles)
dx = (b - a) / n; % Width of each rectangle
area = 0;
for i = 0:n-1
x = a + i * dx; % Left endpoint
height = x^2; % Function value at the left endpoint
area = area + height * dx; % Rectangle area
end
fprintf('The estimated area under y = x^2 is %.4f\n', area);
Let's work with vectors.
This is how you can print all of the elements in a vector:
x = [4:-2:-4];
for k = 1:length(x)
fprintf('Entry %d of vector x is %d \n', k, x(k))
end
Let's use for loops to add two vectors. These vectors are related to the Tower of Hanoi, a mathmetical puzzle found in computer science that is all about moving a stack of disks on a rod onto one of the other rods in order of decreasing size. The rules are that you can only move one disk at a time and that you must start at the top disk. The minimum number of moves required to solve the puzzle is 
, which is what adding these two vectors will result in for 8 disks:
, which is what adding these two vectors will result in for 8 disks:x = (2*(ones(8,1))).^([1:8]')
y = -1*ones(8,1)
s = zeros(length(x),1);
for i = 1:length(x)
s(i,1) = x(i,1) + y(i,1);
end
s
Let's verify that our program calculated this correctly:
x+y
The outputs match.
Let's write a program for calculating dot products. The dot product of two vectors,
and 
, is defined as
Let's code it up:
m = 6;
x = ones(m,1);
y = 2*ones(m,1);
c = 0;
for i = 1:m
c = c + x(i,1)*y(i,1);
end
Let's check our work.
cMAT = dot(x,y) % Built-in function for dot products
c
Let's try a scalar multiplication program. The multiplication of a scalar α and a column vector x is defined as
Let's code it up:
m = 5;
alpha = 16;
x = ones(m,1);
y = zeros(m,1);
yMAT = alpha*x; % Expected output
for i = 1:m
y(i,1) = alpha*x(i,1);
end
Let's check our work. If two things are identical and you subtract one from the other, the answer will of course be zero:
yMAT - y
Our algorithm checks out, but what if we wanted to make it more efficient in terms of memory?
Right now, our code multiplies each entry of x by alpha and then stores it in y. What if we stop keeping track of y and directly write into the existing variable?
m = 5;
alpha = 16;
x = ones(m,1);
yMAT = alpha*x; % Expected output
for i = 1:m
x(i,1) = alpha*x(i,1);
end
Let's check our work again:
yMAT - x
It checks out, with one less variable used. This means that we saved some memory. Using tic and toc, let's see if this makes a difference in time when we set we set m = 1e6:
Indeed, with a much larger number, the time to execute the second version of our algorithm was slightly faster. These differences can add up when dealing with very large programs.
Let's code a SAXPY (single-precision alpha·x plus y) operation:
In other words, 
where x, y, and z are vectors and α is the scalar.
where x, y, and z are vectors and α is the scalar.Let's code it up:
m = 4;
alpha = 8;
x = ones(m,1);
y = randi(10,m,1);
wMAT = alpha*x + y; % Expected output
for i = 1:m
y(i,1) = y(i,1) + alpha*x(i,1);
end
Let's check our work by just putting the vectors next to each other:
[wMAT, y]
If you're familiar with Euler's method, a numerical method for approximating the solution to differential equations, try writing your own program for 
. The solution is below:
. The solution is below:% Define step size h and the number of steps
h = 0.1; % Step size
x0 = 0; % Initial x value
y0 = 1; % Initial y value
x_end = 1; % The final value of x
N = (x_end - x0) / h; % Number of steps
% Initialize arrays to store x and y values
x = zeros(1, N+1);
y = zeros(1, N+1);
x(1) = x0;
y(1) = y0;
for n = 1:N
x(n+1) = x(n) + h; % Compute next x value
y(n+1) = y(n) + h * y(n); % Euler's formula for next y
y(n+1) = round(y(n+1), 4); % Round to 4 decimal places after each step
end
% Display results
disp(' x, y_Euler, Actual e^x, Abs Error');
for n = 1:N+1
actual_value = exp(x(n));
abs_error = abs(actual_value - y(n));
fprintf('%0.4f, %0.4f, %0.4f, %0.4f\n', x(n), y(n), actual_value, abs_error);
end
Let's move on to if statements. This is the general form in pseudocode:
if expression
statements
end
The statements are executed if the expression is true:
num = randi(2)-1;
if num == 1
disp('You got 1')
end
Our value of num makes the condition true. And as you can see, the statement is completed.
Here's another example:
x = [1, 2, 3, 4];
if isscalar(x)
fprintf('x is a scalar')
end
x is of course not a scalar, so it returns logical 0. The statement does not run. It goes right to end. Let's change the condition:
x = [1, 2, 3, 4];
if (isscalar(x) | isvector(x))
fprintf('x is a scalar or vector')
X = [1, 2; 3, 4];
end
Let's test with our new x:
X = [1, 2; 3, 4];
if (isscalar(X) | isvector(X))
fprintf('x is a scalar or vector')
end
The condition is now false, so the statement does not run. Our previous statement printed the message and changed x into a matrix.
An effective application of if statements would be when the program rewards the user for inputting the right value. Here is an example using input:
Entering T into the command window:
Now that we have a basic system for getting the user's input, let's work it into an if statement where I ask a question and the user responds:
Let's run it:
We just created a very small quiz. Great work. You can also do this with buttons (choice):
choice = menu('Select your answer', 'A', 'B', 'C', 'D')
Try writing your own if statement program utilizing buttons.
Let's go back to our quiz. If it were an actual learning exercise, you wouldn't be learning very much if it didn't return anything at all if you were wrong, right? You need some type of feedback. Thankfully, there are ways to address this. First off is the if-else statement. This is the general form in pseudocode:
if condition
statement1
else
statement2
end
If the condition is true, statement1 runs. If it is false, statement2 runs.
Let's adjust our quiz to be more useful:
disp('In the Battle of Britain, the Royal Air Force defeated the Luftwaffe.')
prompt = 'Choose T or F : \n';
user_Choice = input(prompt, 's');
if user_Choice == 'T'
fprintf('Correct!')
else
fprintf('Incorrect. The Luftwaffe being spread too thin and with inferior radio technology led to the first major German defeat of WWII.')
end
Let's run it:
Great, now we have a quiz that actually gives us some feedback.
Here is the last example, which uses nested if statements. It is a "Guess the Number" type game. A nested if statement is when you place one if statement inside another if or else statement:
win = randi(10); % This generates the correct number
menu_prompt = 'Want to play guess the number?'; % This creates the prompt for our menu
user_Choice = menu(menu_prompt, 'Yes', 'No'); % This asks the user if they want to play
if user_Choice == 1 % If the user says yes...
dlg_prompt = 'Guess an integer between 1 and 10.'; %
title = 'Guess the Number';
user_Guess = inputdlg(dlg_prompt, title); % This takes in the user input
guess = str2double(user_Guess); % This converts the string input to a double
if guess == win % This sees if the user guessed the correct number
fprintf('You win the game! \n')
else % This runs if they did not guess the correct number
fprintf('You did not win.')
end
else % If they chose "No" at the start, this statement runs
fprintf('Consider giving it a try next time. Goodbye! \n')
end
Let's run through it:
I will select Yes:
I guessed 5:
I unfortunately did not win, but it looks like our program is working correctly.
Keep in mind that overusing nested if statements can make code hard to read and work with.
Let's move on to elseif statements. This is the general form in pseudocode:
if condition1
statement1
elseif condition2
statement2
elseif condition3
statement3
else
statement4
end
Here's a flowchart:
Try tracing the following example given an input of 65. What clothing will you bring?
if temperature > 90 then
PRINT "It's hot outside. Wear your swimwear."
elseif temperature > 70 then
PRINT "It's warm. Wear a t-shirt."
elseif temperature > 50 then
PRINT "It's chilly. Bring a sweater."
else
PRINT "It's cold. Wear a down jacket."
end
Let's write a program that determines whether the input is a column or row vector:
A = ones(10,1);
B = zeros(1,8);
[mA, nA] = size(A); % Gets number of rows and columns of A
[mB, nB] = size(B); % Gets number of rows and columns of B
% Checks if both inputs are actually vectors
if (~isvector(A)) | (~isvector(B))
fprintf('At least one input is not a vector. \n')
end
% Checks if both inputs are column vectors
if nA == 1 && nB == 1
fprintf('Both inputs are column vectors. \n')
% Checks if both inputs are row vectors
elseif mA == 1 && mB == 1
fprintf('Both inputs are row vectors. \n')
% Runs if they are distinct
else
fprintf('One input is a column vector and the other is a row vector. \n')
end
Let's change A to a row vector to get a different output:
A = ones(1,10); % Changed to a row vector
B = zeros(1,8);
[mA, nA] = size(A);
[mB, nB] = size(B);
% Checks if both inputs are actually vectors
if (~isvector(A)) | (~isvector(B))
fprintf('At least one input is not a vector. \n')
end
% Checks if both inputs are column vectors
if nA == 1 && nB == 1
fprintf('Both inputs are column vectors. \n')
% Checks if both inputs are row vectors
elseif mA == 1 && mB == 1
fprintf('Both inputs are row vectors. \n')
% Runs if they are distinct
else
fprintf('One input is a column vector and the other is a row vector. \n')
end
It looks like our program is running properly. Let's try to improve it by giving the user more specific feedback:
A = ones(10,1);
B = zeros(1,8);
[mA, nA] = size(A);
[mB, nB] = size(B);
% Checks if both inputs are actually vectors
if (~isvector(A)) | (~isvector(B))
fprintf('At least one input is not a vector. \n')
end
% Checks if both inputs are column vectors
if nA == 1 && nB == 1
fprintf('A is a column vector and B is a column vector. \n')
% Checks if both inputs are row vectors
elseif nA == 1 && mB == 1
fprintf('A is a column vector and B is a row vector. \n')
% Runs if they are distinct
elseif mA == 1 && nB == 1
fprintf('A is a row vector and B is a column vector. \n')
else
fprintf('A is a row vector and B is a row vector. \n')
end
Let's explain the concept of short-circuit logical operators. Recall the AND and OR operators.
For AND, let's say have two parameters. The first input is false. The entire expression will always be false, no matter what the second input is.
For OR, let's say we have two parameters. The first input is true. The entire expression will always be true, no matter what the second input is.
Thus, there is no point in checking the second input here. Short-circuit AND (&&) and Short-circuit OR (||) will stop evaluating once the result is determined. This means more efficient code. Here is a pseudocode example of where this could be helpful.
This is the code block that is ran after the user logs in:
user_is_logged_in = true;
if user_is_logged_in | verify_login
grant access
end
There's no point in testing verify_login. The overall statement will be true regardless of its result. Why would a user who is already logged in need to verify themselves again?
user_is_logged_in = true;
if user_is_logged_in || verify_login
grant access
end
Now with the short-circuit OR (||), the verify_login check will be skipped to save resources.
9. Writing Programs Using Function Files
You can add custom functions to MATLAB using .m files. This is the general syntax:
function(output1, output2, ..., outputn) = myFunction(input1, input2, inputm)
The number of inputs does not need to equal the number of outputs. You may wonder how function files differ from script files as they both use the .m extension. Both allow you to reuse sequences of code. Scripts are much more simple, as they store commands in the same way you would enter them into the command window. However, functions are more versatile as they have their own environment. They can pass input values and return output values. They are for extending your library rather than just automating a series of pre-existing functions.
Remember that when creating something like this, you want it to be reusable, scaleable, and as convenient as possible for your peers. Here are some practices and style rules you want to follow:
- Code should have formal documentation: what the code does, how it works, the functions it uses, any other files it requires, a quickstart guide (like an example), and your information.
- Write code that is easy for humans to read and not just computers. Don't write code that is only good enough to be ran once, published in Cell Press, and then linked as a footnote that no one will ever use again.
- The names you use for variables or functions should be consistent (like using camelCase), descriptive, and unique.
- For these variables, make sure to use relevant and consistent data types.
- Make it easy for people to share and cite your work.
- It's generally best to use established and open-source libraries.
- Minimize redundant code
- Indentation should be consistent and used to make loops discernible.
- Use commas in arrays and follow them up with spaces.
- Surround operators with spaces
- Do not submit code with commented out functions or TODOs.
I've linked these guides in "Further Reading." Some are not necessarily for MATLAB, but the ideas still apply. If your professor or group has specific rules, you should obviously follow what they want. I apologize for any errors or inconsistencies in this guide.
Now, we will start working with function files. Throughout this section, you will learn a little bit about program development and linear algebra.
Let's create our first function that determines if two vectors are the same size. First, let's make a new directory for it:
mkdir vectorSizeTest
cd vectorSizeTest
Next, go to New and select Function. We now have our file. This is what MATLAB has provided us by default to start:
function [outputArg1,outputArg2] = untitled(inputArg1,inputArg2)
%UNTITLED Summary of this function goes here
% Detailed explanation goes here
outputArg1 = inputArg1;
outputArg2 = inputArg2;
end
Let's also create a test driver file called vectorSizeTestDriver.m.
Here is the function file:
function [vectorTest] = vectorSizeTest(x,y)
% This function tests to see if two input vectors have the same number of elements
%
% Inputs:
% x — input vector
% y — input vector
% where x and y are nonempty vectors
%
% Outputs:
% vectorTest — logical scalar to tell if inputs have same number of elements
%
% Checks that user is passing correct number of inputs
if nargin ~= 2
error('* This function must have two inputs. *')
return
end
% Checks that neither input is empty
if isempty(x) || isempty(y)
error('* At least one input is empty. *')
return
end
% Checks that both inputs are vectors
if (~isvector(x)) || (~isvector(y))
error('* At least one input is not a vector. *')
end
% Get the dimensions of input x
[mx, nx] = size(x);
% Get the dimensions of input y
[my, ny] = size(y);
if (mx == my) && (nx == ny) % If both are same sized column or row vectors
vectorTest = true;
return
elseif (mx == ny) && (nx == my)
vectorTest = true;
return
else
vectorTest = false;
return
end
Let's do some manual test driving in our script file.
For empty arrays, this is the error we created:
This tests to make sure we have same-sized arrays:
Do not include clear,clc in any code you send to someone else if it is going to delete all their progress. It's only there for us as we are in the design process.
Let's make another function file for scalar multiplication. Remember that we did this before for a script file, but it was very rudimentary.
The multiplication of a scalar α and a column vector x is defined as:
The pseudocode for this:
for i = 1:m
y(i,1) = alpha*x(i,1)
end
The multiplication of a scalar α and a row vector x is defined as:
The pseudocode for this:
for k = 1:m
y(1,k) = alpha*x(1,k)
end
Before we start coding, let's think about what we want to do:
- Takes in two inputs: a scalar α and a vector x stored in double precision
- If the user does not supply the proper number of inputs, an error message is produced.
- If either inputs are empty, an error message is produced.
- If scalar α is not a scalar or if vector x is not a vector, an error message is produced.
- It will output the product
. - The input and output should have the same dimensions.
Let's get started by creating a folder called BLAS1_Library for basic linear algebra subprograms level 1. Then, add two files to to it called dscal_test_driver.m and dscal.m. For this program, we are once again kind of working backwards: starting with error testing and then moving into the main function. This isn't necessarily the best practice, but it works for this purpose.
Here is our current dscal.m:
function xOut = dscal(alpha,x)
% dscal scales a vector x by a scalar alpha.
%
% Vector x can be either a column or row vector.
% The output xOut is a vector with dimensions identical to input vector x.
% All data should be stored in double precision.
%
% Syntax:
% xOut = dscal(alpha,x) = alpha*x
%
% Inputs:
% alpha — scalar alpha
% x — array x
%
% Outputs:
% xOut
%
%
% Any other required files, tools, or software:
%
% See also/references:
%
% Author:
% Work/school email:
% Website:
% Created on:
% Last revision:
% Version:
% Copyright information:
% Checks that user is passing correct number of inputs
if nargin ~= 2
error('* This function must have two inputs. * ')
end
% Checks that neither input is empty
if isempty(alpha) || isempty(x)
error('* At least one input is empty. *')
end
% Checks that both inputs are stored in double precision
if ~isa(alpha, 'double') || ~isa(x, 'double')
error('* Both inputs must be stored in double precision. *')
end
% Checks that first input is stored as a scalar
if ~isscalar(alpha)
error('* Input alpha must be a scalar. *')
end
% Gets the dimensions of input x
[mx,nx] = size(x)
% Checks that second input is a vector
if mx ~= 1 && nx ~= 1
error('* Second input x must be a vector. *')
end
And here is our test driving in dscal_test_driver.m:
% This script file is used to test drive the dscal.m function
% ------------------ TEST DRIVING ------------------
% % Error checking
% If user doesn't provide exactly two inputs
% alpha = 5;
% x = ones(15,1);
% [xOut] = dscal(alpha);
%
% If user provides an empty input
% alpha = [];
% x = ones(15,1);
% [xOut] = dscal(alpha,x);
%
% If user doesn't use the double precision data type
% alpha = double(-3);
% x = uint8(ones(15,1));
% [xOut] = dscal(alpha,x);
%
% Checks that first input alpha is a scalar
% alpha = [1,2];
% x = ones(15,1);
% [xOut] = dscal(alpha,x);
%
% Checks that second input x is a vector
% alpha = 2;
% x = ones(2,14);
% [xOut] = dscal(alpha,x);
I went one by one for each of these cases and then commented them out when I was done.
Our error testing is done. Let's move into the main algorithm. Here is the rest of dscal.m:
% Main algorithm
if nx == 1 % For column vectors
for i = 1:mx
x(i,1) = alpha*x(i,1);
end
else
for i = 1:nx
x(1,i) = alpha*x(1,i);
end
end
xOut = x;
return
end
Here is the corresponding test driving:
% % Verify Output
% Multiply scalar alpha times a column vector
m = 6;
alpha = randi(16, 1, 1);
x = randi(21,m,1) - 10*ones(m, 1);
xOut = dscal(alpha,x);
fprintf('Comparison of MATLAB output to our function output. \n')
[alpha*x, xOut]
% Multiply scalar alpha times a row vector
n = 10;
alpha = randi(16, 1, 1);
x = randi(21, 1, n) - 10*ones(1, n);
xOut = dscal(alpha,x);
fprintf('Comparison of MATLAB output to our function output. \n')
[alpha*x; xOut]
It looks like our algorithm is working properly:
Let's make another program that copies a vector into a different vector. Let's start by duplicating our previous files and naming the copies dcopy_test_driver.m and dcopy.m.
Here is our complete test driving:
% This script file is used to test drive the dcopy.m function
% ------------------ TEST DRIVING ------------------
% % Error checking
% If user doesn't provide exactly two inputs
% x = ones(15,1);
% y = ones(15,1);
% [yOut] = dcopy(x);
%
% If user provides an empty input
% x = ones(15,0);
% y = ones(15,1);
% [yOut] = dcopy(x,y);
%
% If user doesn't use the double precision data type
% x = single(ones(15,1));
% y = single(ones(15,1));
% [yOut] = dcopy(x,y);
%
% Checks that second input x is a vector
% x = ones(15,2);
% y = ones(1,10);
% [yOut] = dcopy(x,y);
%
% Checks that both inputs have the same number of entries
% x = ones(10,1);
% y = ones(4,1);
% [yOut] = dcopy(x,y);
% % Verify Output
% Copy a column vector onto a column vector
m = 8;
x = randi(21,m,1) - 10*ones(m, 1);
y = randi(21,m,1) - 10*ones(m, 1);
yOut = dcopy(x,y);
fprintf('Comparison of original column vector with copy column vector \n')
[x, yOut]
% Copy a row vector onto a row vector
m = 8;
x = randi(21,1,m) - 10*ones(1, m);
y = randi(21,1,m) - 10*ones(1, m);
yOut = dcopy(x,y);
fprintf('Comparison of original row vector with copy row vector \n')
[x; yOut]
% Copy a column vector onto a row vector
m = 8;
x = randi(21, m, 1) - 10*ones(m, 1);
y = randi(21, 1, m) - 10*ones(1, m);
yOut = dcopy(x, y);
fprintf('Comparison of original column vector with copy row vector \n')
disp(x)
disp(yOut)
% 4. Row → Column
m = 8;
x = randi(21, 1, m) - 10*ones(1, m);
y = randi(21, m, 1) - 10*ones(m, 1);
yOut = dcopy(x, y);
fprintf('Comparison of original row vector with copy column vector \n')
disp(x)
disp(yOut)
Here is our function file:
function yOut = dcopy(x,y)
% dcopy copies a vector x into a vector y.
%
% Vector x can be either a column or row vector.
% All data should be stored in double precision.
% Inputs x and y should have the same number of entries and be nonempty.
% The dimensions of xOut should be identical to those of y
% A row vector input should output a row vector.
% A column vector input should output a column vector.
%
% Syntax:
% yOut = dcopy(x,y)
%
% Inputs:
% x - vector that will be scaled
% y — vector that determines size of output
%
% Outputs:
% yOut
%
% Any other required files, tools, or software:
%
% See also/references:
%
% Author:
% Work/school email:
% Website:
% Created on:
% Last revision:
% Version:
% Copyright information:
% Checks that user is passing correct number of inputs
if nargin ~= 2
error('* This function must have two inputs. * ')
end
% Checks that neither input is empty
if isempty(x) || isempty(y)
error('* At least one input is empty. *')
end
% Checks that both inputs are stored in double precision
if ~isa(x, 'double') || ~isa(y, 'double')
error('* Both inputs must be stored in double precision. *')
end
% Gets the dimensions of input x
[mx,nx] = size(x);
% Gets the dimensions of input y
[my,ny] = size(y);
% Checks that both inputs are vector
if (mx ~= 1 && nx ~= 1) || (my ~= 1 && ny ~= 1)
error('* Both inputs must be vectors. *')
end
% Checks that inputs x and y have the same length
if (mx*nx ~= my*ny)
error('* Both inputs must have the same number of elements. *')
end
% Main algorithm
if nx == 1 % x is a column vector
if ny == 1 % y is a column vector
for i = 1:mx
y(i,1) = x(i,1);
end
else % y is a row vector
for i = 1:mx
y(1,i) = x(i,1);
end
end
else % x is a row vector
if ny == 1 % y is a column vector
for i = 1:nx
y(i,1) = x(1,i);
end
else % y is a row vector
for i = 1:nx
y(1,i) = x(1,i);
end
end
end
yOut = y;
return
end
Notice how we were able to use our previous test driver and function file as templates for our current one.
Let's test it out:
Awesome!
Let's go through a more formal outline for program development.
- Understand the problem: learn any relevant context and analyze the problem
- Come up with a general solution: design algorithms and create program specifications
- Implement: convert your algorithm into code, test, and debug
- Upkeep: fix any bugs that are discovered, add new features, and optimize
- Grow: repurpose your program for other applications, distribute, and acquire a license if you desire (hopefully an open source one!)
Let's write a function for multiplying a matrix by a column vector. Let's follow step 1 and make sure we understand the problem. First, the definition:
Let's write that out using colon notation:
Alternatively, we can write it as a summation:
A linear combination is an expression formed by adding up a set of vectors where each is multiplied by a constant.
Here is an example using linear combinations:
, 
Using our definition,
Keep this in mind:
The inner dimensions must match up.
More commonly, this whole operation can be can be seen as dot products. Every column of A gets multiplied by a vector x to produce a new vector b. Thus, the number of columns in A must match the number of entries given by x. The vector b will be in this form:
where each element is a dot product
For instance, 
is the dot product of 
and 
. We will then loop through each column of A, multiplying it by the corresponding element in x to form b.
is the dot product of and . We will then loop through each column of A, multiplying it by the corresponding element in x to form b.Here it is again:
This yields the same result:
Sometimes, if you are working on a problem in one language, it's a good idea to try it in another. Let's verify our example using Python:
import numpy as np
A = np.array([
[1, 2, 3],
[4, 5, 6],
[7, 8, 9],
[10, 11, 12]
])
x = np.array([
[-5],
[-10],
[-15]
])
b = A @ x
print(b)
It outputs this:
[[-70]
[-160]
[-250]
[-340]]
This matches our above work.
Let's repeat but this time for row vectors:
Here it is in colon notation and as a summation:
Let's look at an example:
, 
Here's the calculation:
Now, let's do it with dot products:
Here is another visual representation:
= 
Notice how that was a fairly elaborate explanation of an easy calculation. It is common for engineers to spend much more time thinking about problems than actually building.
Now that step 1 is done, let's move on to step 2. We understand the problem, so now it's time to come up with a solution. Let's think about what information is relevant to solving it.
We have two types of multiplication:
Matrix-column-vector multiplication
- Vector version, where data are chunked into vector-sized pieces
- Scalar version, where data are cut up into individual entries
Row-vector-matrix multiplication
- Vector version, where data are chunked into vector-sized pieces
- Scalar version, where data are cut up into individual entries
Vector Version | Scalar version | |
Matrix-column-vector multiplication | Column partition   | Cut up into individual pieces  |
Row vector matrix multiplication | Row partition   | Cut up into individual pieces  |
How do we distinguish between these options? When do we want them? Do they matter?
Partitioning matrices has many applications, like Strassen's algorithm, parallelism, and ray tracing. They can let us divide up and analyze data in more efficient ways.
It's all about understanding the different angles available for a problem.
Let's take a look at some of the propositions we have so far:
Proposition 1 for matrix-column-vector multiplication:
Multiplying a matrix by a column vector and then taking the transposing its result is the same as transposing the column vector and matrix.
Let's prove it:
We have
Taking the transpose of
Now for 
Since scalar multiplication is commutative, the two expressions are equal:
Let's prove it in reverse:
Take the transpose of both sides:
By the property of double transposes, we have this:
Thus, proposition 1 holds. This is because taking the transpose twice returns the original matrix.
Proposition 2 for row-vector-matrix multiplication:
Multiplying a row vector by a matrix and then transposing the result is the same as multiplying the matrix by the original column vector.
Let's prove it:
We have 
Compute the 
Now compute the 
Since scalar multiplication is commutative, the two expressions are equal:
Let's prove it in reverse:
Take the transpose of both sides:
By the property of double transposes, we have this:
Thus, proposition 2 holds.
It's important that we spent so much time going over the different ways to approach this problem because now we have something useful we can apply to our coding:
Let's say we have an instance row-vector-matrix multiplication:
Knowing that row-vector-matrix multiplication is equivalent to the transpose of column vector-matrix multiplication, we can do the latter first. This is helpful if say, we didn't want to deal with row vectors:
This gives us a column vector which is the transpose of the result we want.
Thus, we have this:
We computed a row-vector-matrix product using only column vectors. In other words, we did it using a single function.
It's often better to reuse or repurpose a function than writing brand new functions every time. This ties into the idea of modular programming.
Let's write an algorithm for matrix-column-vector multiplication using linear combinations:
y = zeros(m,1);
for k = 1:n
y = y + x(k,1) * A(:,k)
end
Let's implement it into MATLAB using a previous example. Notice how this isn't our final function. It's what's known as a playground or sandbox, which is a local environment where you can test your code:
A = [1, 2, 3;
4, 5, 6;
7, 8, 9;
10, 11, 12];
x = [-5; -10; -15];
yMAT = A*x;
[m,n] = size(A);
y = zeros(m,1);
for k = 1:n
y = y + x(k,1) * A(:,k);
end
y
Here's an an alternative way of looking at this:
Here's iteration k = 1:
x_1 = -5
A(:,1) = [1; 4; 7; 10]
x_1 * A(:,1)
= [-5;
-20;
-35;
-50]
y = y + x_1 * A(:,1)
= [0;
0;
0;
0]
+ [-5;
-20;
-35;
-50]
= [-5;
-20;
-35;
-50]
Here's iteration k = 2:
x_2 = -10
A(:,2) = [2; 5; 8; 11]
x_2 * A(:,2)
= [-20;
-50;
-80;
-110]
y = y + x_2 * A(:,2)
= [-5;
-20;
-35;
-50]
+ [-20;
-50;
-80;
-110]
= [-25;
-70;
-115;
-160]
Here's iteration k = 3:
x_3 = -15
A(:,3) = [3; 6; 9; 12]
x_3 * A(:,3)
= [-45;
-90;
-135;
-180]
y = y + x_3 * A(:,3)
= [-25;
-70;
-115;
-160]
+ [-45;
-90;
-135;
-180]
= [-70;
-160;
-250;
-340]
This looks a bit like a SAXPY operation, which was previously used as an example. Let's implement this method:
A = [1, 2, 3;
4, 5, 6;
7, 8, 9;
10, 11, 12];
x = [-5; -10; -15];
yMAT = A*x;
[m,n] = size(A);
y = zeros(m,1);
for k = 1:n
for i = 1:m
y(i,1) = y(i,1) + x(k,1) * A(i,k);
end
end
y
This algorithm iterates entry-by-entry instead. It still accomplishes the same task.
Let's do it using dot products formally—back to the drawing board!
Here's the algorithm:
y = zeros(m,1);
for i = 1:m
y(i,1) = y(i,1) + dot(A(i,:), x)
end
Let's implement it:
A = [1, 2, 3;
4, 5, 6;
7, 8, 9;
10, 11, 12];
x = [-5; -10; -15];
yMAT = A*x;
[m,n] = size(A);
y = zeros(m,1);
for i = 1:m
y(i,1) = y(i,1) + dot(A(i,:),x);
end
y
Let's do this same process but for row-vector-matrix multiplication. First, let's write the algorithm:
y = zeros(1,n);
for i = 1:m
y = y + x(i,1) * A(i,:)
end
Let's implement it:
A = [1, 5, 9;
2, 6, 10;
3, 7, 11;
4, 8, 12];
x = [-2; -4; -6; -8];
yMAT = x'*A;
[m,n] = size(A);
y = zeros(1,n);
for i = 1:m
y = y + x(i,1) * A(i,:);
end
y'
We can do it again but iterating entry-by-entry instead:
A = [1, 5, 9;
2, 6, 10;
3, 7, 11;
4, 8, 12];
x = [-2; -4; -6; -8];
yMAT = x'*A;
[m,n] = size(A);
y = zeros(1,n);
for i = 1:m
for k = 1:n
y(1,k) = y(1,k) + x(i,1) * A(i,k);
end
end
y'
This type of testing that we are doing is pretty cursory. We can test using random and much larger matrices too:
A = rand(50,5);
x = rand(50,1);
yMAT = x'*A;
[m,n] = size(A);
y = zeros(1,n);
for i = 1:m
for k = 1:n
y(1,k) = y(1,k) + x(i,1) * A(i,k);
end
end
error = norm(yMAT - y)
Our error is a very small number due to how MATLAB rounds integers. It is essentially zero.
Let's repeat the process using dot products. Here is the algorithm:
y = zeros(1,n);
for k = 1:n
y(1,k) = y(1,k) + dot(x, A(:,k))
end
Let's implement it:
A = [1, 5, 9;
2, 6, 10;
3, 7, 11;
4, 8, 12];
x = [-2; -4; -6; -8];
yMAT = x'*A;
[m,n] = size(A);
y = zeros(1,n);
for k = 1:n
y(1,k) = y(1,k) + dot(x, A(:,k));
end
y'
Before we produce our function, let's summarize the data flow for our algorithms. We can do this using tables.
Here's matrix-column-vector multiplication:
Loop order | Operation in outer loop | Inner loop data access | Purpose | Algorithm name |
K i | SAXPY | A by columns, x by rows | Adds a scaled column of A to y | column-SAXPY |
i k | dot product | A by rows, x by rows | Computes a scalar entry of y by dotting a row of A with x | row-dot |
Here's row vector-matrix multiplication:
Loop order | Operation in outer loop | Inner loop data access | Purpose | Algorithm name |
i k | SAXPY | A by columns,  by rows | Adds a scaled row of A to y | row-SAXPY |
k i | dot product | A by rows, by columns | Computes a scalar entry by dotting a column of A with x | column-dot |
Let's add to why thinking about this is important:
Matrix-vector algorithms interact with a computer's memory. Most CPUs follow the fetch-decode-execute cycle, an idea from John von Neumann.
- Fetch the instruction from memory
- Decode the instruction and any necessary data into a language the CPU understands
- Execute and save if tasked
This process then repeats. Thus, loop order and the ways a program accesses memory can affect how efficiently the CPU fetches data. You may not see much of a difference in your work for a while, but more advanced programs can be hundreds to thousands of times larger than yours. These types of things make a difference over time.
Let's start making some specifications for our function.
First, here's the function:
where:
This is a subroutine from BLAS Level 2 that was written in 1986. Indeed, we are recreating part of LAPACK, the library upon which MATLAB itself was built.
Now, let's write out the specs. Our function should:
Take in Trans, which indicates if we want to multiply using 
- If Trans = 'N', then the output should be
. - If Trans = 'T', then the output should be
.
Take in two scalars α and β stored in double precision
- If input beta is set to 0, then y is irrelevant and can be set to 0.
- If input beta is nonzero, then input y must be provided.
Take as input an matrix A stored in double precision
matrix A stored in double precisionTake as input two column vectors x and y stored in double precision
- If Trans = 'N', then x is
and y is . - If Trans = 'T', then x is and y is .
An error message is produced if:
- The user does not use the proper number of inputs.
- Trans ≠ 'N' or 'T'
- If either α or β is not a scalar stored in double precision
- If A is not a matrix stored in double precision
- If either x or y are not vectors with proper dimension stored in double precision
Make a new folder BLAS2_Library. Create a script file dgemv_test_driver.m and function file dgemv.m. Reuse the documentation from our BLAS1_Library.
Here's our documentation and the first error check in dgemv.m:
function yOut = dgemv(Trans, alpha, A, x, beta, y)
% dgemv performs a matrix-vector multiplication in the form:
% y := alpha*A*x + beta*y or
% y := alpha*A^T*x + beta*y
% where alpha and beta are scalars,
% x and y are column vectors, and
% A is an m-by-n matrix.
%
% Syntax:
% yOut = dgemv('N', alpha, x, A, beta, y)
% yout = dgemv('TT', alpha, x, A, beta, y)
%
% Inputs:
% Trans: stored as character data
% specifies if we want to calculate using the transpose
% Trans = 'N' (no transpose) means y := alpha*A*x + beta*y
% Trans = 'T' (transpose) means y := alpha*A^t*x + beta*y
% This variable is case sensitive.
%
%
% alpha: scalar alpha stored in double precision
%
% A: two dimensional array A stored in double precision
% stores entries of matrix used for the matrix-vector product
% has dimensions m-by-n
%
% x: one dimensional array x stored in double precision
% stores entries of column vector x
% when trans = 'N', x must have dimensions n-by-1.
% when trans = 'T', x must have dimensions m-by-1.
%
% beta: scalar beta stored in double precision
% if beta is set to zero, then y does not need to be defined.
%
% y: one dimensional array y stored in double precision
% if input beta is nonzeo, array stores column vector y.
% when Trans = 'N', y must have dimensions m-by-1.
% when trans = 'T', y must have dimensions n-by-1.
% y is overwritten by updated vector y.
%
% Outputs:
% yOut: if Trans = 'N', then yOut := alpha*A*x + beta*y.
% if Trans = 'T', then yOut := alpha*A^T*x + beta*y.
% stored in double precision
% this is the desired matrix-vector product.
%
% Any other required files, tools, or software: dgemv_test_driver.m
%
% See also/references:
%
% Author:
% Work/school email:
% Website:
% Created on:
% Last revision:
% Version:
% Copyright information:
% Check that user passes five inputs into the function
if nargin <5
error('* You must have at least 5 inputs. *')
end
yOut = y;
return
Let's run it in dgemv_test_driver.m:
% % Error checking
% If user doesn't provide five inputs
Trans = 'N';
alpha = 1;
A = [1, 2, 3; 4, 5, 6; 7, 8, 9; 10, 11, 12];
x = [-1; -2; -3];
yOut = dgemv(Trans, alpha, A, x);
We only have 4 inputs, so it returns errors as expected:
Similarly to the process in previous function files, we can now comment out this error checking instance in the test driver. You can comment out and uncomment multiple lines of code at once in the Editor tab.
Keep on repeating this process for the rest of the specs.
Here is the rest of the error checking statements in dgemv.m:
% Checks that the first input is valid
if Trans ~= 'N' && Trans ~= 'T'
error('* First input must be a character equal to N or T. Input is case sensitive. *')
end
% Checks that alpha and beta are nonempty scalars stored in double precision
if ~isscalar(alpha) || ~isscalar(beta)
error('* Inputs alpha and beta must be stored as nonempty scalars. *')
elseif ~isa(alpha, 'double') || ~isa(beta, 'double')
error(' Inputs alpha and beta must be stored as doubles. *')
end
% Checks that input A is a nonempty array stored in double precision
if isempty(A)
error('* Input array A must be nonempty. *')
elseif ~ismatrix(A)
error('* Input array A must be a matrix. *')
elseif ~isa(A, 'double')
error(' Input array A must be stored as a double. *')
end
% Checks handling of beta and vector y
% Get dimensions of matrix A
[m,n] = size(A);
% Set lenx and leny to be the lengths of vectors x and y
if Trans == 'N'
lenx = n;
leny = m;
elseif Trans == 'T'
lenx = m;
leny = n;
end
% Form vector y = beta*y
if beta == 0
y = zeros(leny, 1);
elseif beta ~= 0 && nargin ~= 6
error(' * When beta is nonzero, vector y must be given as input. *')
elseif beta ~= 0 && any(size(y) ~= [leny, 1])
error(' * Input vector y does not have proper dimensions. *')
else
for i = 1:leny
y(i,1) = beta*y(i,1);
y
end
end
% If alpha = 0, return yOut = beta*y;
if alpha == 0
yOut = y;
return
end
% Checks input x is nonempty and stored in double precision
if any(size(x) ~= [lenx, 1])
error('* Input vector x does not have proper dimensions. *')
elseif ~isa(x, 'double')
error('* Input x must be stored as a double. *')
end
Here is the rest of the error test driving in dgemv_test_driver.m:
% Check that the first input is 'N' or 'T'
% Trans = 'C';
% alpha = 1;
% A = [1, 2, 3; 4, 5, 6; 7, 8, 9; 10, 11, 12];
% x = [-1; -2; -3];
% beta = 0;
% yOut = dgemv(Trans, alpha, A, x, beta);
%
% Check that alpha and beta are nonempty scalars in double precision
% Trans = 'N';
% alpha = [];
% A = [1, 2, 3; 4, 5, 6; 7, 8, 9; 10, 11, 12];
% x = [-1; -2; -3];
% beta = double(0);
% yOut = dgemv(Trans, alpha, A, x, beta);
%
% Check that A is a nonempty array in double precision
% Trans = 'N';
% alpha = 1;
% A = zeros(10, 10, 3);
% x = [-1; -2; -3];
% beta = 0;
% yOut = dgemv(Trans, alpha, A, x, beta);
%
% Check handling of nonzero input scalar beta and vector y
% Trans = 'N';
% alpha = 1;
% A = [1, 2, 3; 4, 5, 6; 7, 8, 9; 10, 11, 12];
% x = [-1; -2; -3];
% beta = 1;
% yOut = dgemv(Trans, alpha, A, x, beta);
%
% Check for errors in input vector x
% Trans = 'N';
% alpha = 1;
% A = [1, 2, 3; 4, 5, 6; 7, 8, 9; 10, 11, 12];
% x = single([-1; -2; -3]);
% beta = 2;
% y = ones(4,1);
% yOut = dgemv(Trans, alpha, A, x, beta, y);
We are now ready to start coding the main algorithm.
After coding it up, here is our final dgemv.m:
function yOut = dgemv(Trans, alpha, A, x, beta, y)
% dgemv performs a matrix-vector multiplication in the form:
% y := alpha*A*x + beta*y or
% y := alpha*A^T*x + beta*y
% where alpha and beta are scalars,
% x and y are column vectors, and
% A is an m-by-n matrix.
%
% Syntax:
% yOut = dgemv('N', alpha, x, A, beta, y)
% yout = dgemv('TT', alpha, x, A, beta, y)
%
% Inputs:
% Trans: stored as character data
% specifies if we want to calculate using the transpose
% Trans = 'N' (no transpose) means y := alpha*A*x + beta*y
% Trans = 'T' (transpose) means y := alpha*A^t*x + beta*y
% This variable is case sensitive.
%
%
% alpha: scalar alpha stored in double precision
%
% A: two dimensional array A stored in double precision
% stores entries of matrix used for the matrix-vector product
% has dimensions m-by-n
%
% x: one dimensional array x stored in double precision
% stores entries of column vector x
% when trans = 'N', x must have dimensions n-by-1.
% when trans = 'T', x must have dimensions m-by-1.
%
% beta: scalar beta stored in double precision
% if beta is set to zero, then y does not need to be defined.
%
% y: one dimensional array y stored in double precision
% if input beta is nonzeo, array stores column vector y.
% when Trans = 'N', y must have dimensions m-by-1.
% when trans = 'T', y must have dimensions n-by-1.
% y is overwritten by updated vector y.
%
% Outputs:
% yOut: if Trans = 'N', then yOut := alpha*A*x + beta*y.
% if Trans = 'T', then yOut := alpha*A^T*x + beta*y.
% stored in double precision
% this is the desired matrix-vector product.
%
% Any other required files, tools, or software: dgemv_test_driver.m
%
% See also/references:
%
% Author:
% Work/school email:
% Website:
% Created on:
% Last revision:
% Version:
% Copyright information:
% Checks that user passes five inputs into the function
if nargin <5
error('* You must have at least 5 inputs. *')
end
% Checks that the first input is valid
if Trans ~= 'N' && Trans ~= 'T'
error('* First input must be a character equal to N or T. Input is case sensitive. *')
end
% Checks that alpha and beta are nonempty scalars stored in double precision
if ~isscalar(alpha) || ~isscalar(beta)
error('* Inputs alpha and beta must be stored as nonempty scalars. *')
elseif ~isa(alpha, 'double') || ~isa(beta, 'double')
error(' Inputs alpha and beta must be stored as doubles. *')
end
% Checks that input A is a nonempty array stored in double precision
if isempty(A)
error('* Input array A must be nonempty. *')
elseif ~ismatrix(A)
error('* Input array A must be a matrix. *')
elseif ~isa(A, 'double')
error(' Input array A must be stored as a double. *')
end
% Checks handling of beta and vector y
% Get dimensions of matrix A
[m,n] = size(A);
% Set lenx and leny to be the lengths of vectors x and y
if Trans == 'N'
lenx = n;
leny = m;
elseif Trans == 'T'
lenx = m;
leny = n;
end
% Form vector y = beta*y
if beta == 0
y = zeros(leny, 1);
elseif beta ~= 0 && nargin ~= 6
error(' * When beta is nonzero, vector y must be given as input. *')
elseif beta ~= 0 && any(size(y) ~= [leny, 1])
error(' * Input vector y does not have proper dimensions. *')
else
for i = 1:leny
y(i,1) = beta*y(i,1);
y
end
end
% If alpha = 0, return yOut = beta*y;
if alpha == 0
yOut = y;
return
end
% Checks input x is nonempty and stored in double precision
if any(size(x) ~= [lenx, 1])
error('* Input vector x does not have proper dimensions. *')
elseif ~isa(x, 'double')
error('* Input x must be stored as a double. *')
end
% Main algorithm
% Do the operation where the elements of A are accessed through column
% partition in one pass.
if Trans == 'N' % In the form y = alpha*A*x + y
% Matrix-column-vector multiplication via linear combinations
for k = 1:n
temp = alpha*x(k,1);
for i = 1:m
y(i,1) = y(i,1) + temp*A(i,k);
end
end
else % In the form y = alpha*A^T*x + y
% Matrix-column-vector multiplication via dot products
for k = 1:n
temp = 0;
for i = 1:m
temp = temp + A(i, k)*x(i,1);
end
y(k,1) = y(k,1) + temp;
end
end
yOut = y;
return
Here is our final dgemv_test_driver.m, now updated with output examples:
% This script file is used to test drive the dgemv.m function
%
% An error message is produced if:
% The user does not use the proper number of inputs.
% Trans 'N' or 'T'
% If either or is not a scalar stored in double precision
% If is not a matrix stored in double precision
% If either or are not vectors with proper dimension stored in double precision
% ------------------ TEST DRIVING ------------------
% % Error checking
% If user doesn't provide at least five inputs
% Trans = 'N';
% alpha = 1;
% A = [1, 2, 3; 4, 5, 6; 7, 8, 9; 10, 11, 12];
% x = [-1; -2; -3];
% yOut = dgemv(Trans, alpha, A, x);
%
% Check that the first input is 'N' or 'T'
% Trans = 'C';
% alpha = 1;
% A = [1, 2, 3; 4, 5, 6; 7, 8, 9; 10, 11, 12];
% x = [-1; -2; -3];
% beta = 0;
% yOut = dgemv(Trans, alpha, A, x, beta);
%
% Check that alpha and beta are nonempty scalars in double precision
% Trans = 'N';
% alpha = [];
% A = [1, 2, 3; 4, 5, 6; 7, 8, 9; 10, 11, 12];
% x = [-1; -2; -3];
% beta = double(0);
% yOut = dgemv(Trans, alpha, A, x, beta);
%
% Check that A is a nonempty array in double precision
% Trans = 'N';
% alpha = 1;
% A = zeros(10, 10, 3);
% x = [-1; -2; -3];
% beta = 0;
% yOut = dgemv(Trans, alpha, A, x, beta);
%
% Check handling of nonzero input scalar beta and vector y
% Trans = 'N';
% alpha = 1;
% A = [1, 2, 3; 4, 5, 6; 7, 8, 9; 10, 11, 12];
% x = [-1; -2; -3];
% beta = 1;
% yOut = dgemv(Trans, alpha, A, x, beta);
%
% Check for errors in input vector x
% Trans = 'N';
% alpha = 1;
% A = [1, 2, 3; 4, 5, 6; 7, 8, 9; 10, 11, 12];
% x = single([-1; -2; -3]);
% beta = 2;
% y = ones(4,1);
% yOut = dgemv(Trans, alpha, A, x, beta, y);
% % Verify Output
% Matrix-column-vector multiplication
% Multiply matrix A by column vector x with yOut = A*x
% Trans = 'N';
% alpha = 1;
% A = [1, 2, 3; 4, 5, 6; 7, 8, 9; 10, 11, 12];
% x = double([-1; -2; -3]);
% beta = 0;
% yOut = dgemv(Trans, alpha, A, x, beta);
% yMAT = A*x;
% [yMAT, yOut]
% err = yMAT - yOut
% % No error! Our function must be identical to what MATLAB's code does.
% Multiply matrix A by column vector x with yOujt = alpha*A*x + beta*y
% Trans = 'N';
% alpha = 1;
% A = [1, 2, 3; 4, 5, 6; 7, 8, 9; 10, 11, 12];
% x = double([-1; -2; -3]);
% beta = 3;
% y = [1; -1; 1; -1];
% yOut = dgemv(Trans, alpha, A, x, beta, y);
% yMAT = alpha*A*x + beta*y;
% [yMAT, yOut]
% err = norm(yMAT - yOut)
% Row-vector-matrix multiplication
% Multiply matrix A by row vector x^T with yOut = A^T*x
% Trans = 'T';
% alpha = 1;
% A = [1, 2, 3; 4, 5, 6; 7, 8, 9; 10, 11, 12];
% x = double([-1; -2; -3; -4]);
% If we replace A and x with large random matrices, we get very small
% deviation between our function and MATLAB. This is from the rounding
% due to different order of operations or statements.
% beta = 0;
% yMAT = A'*x;
% yOut = dgemv(Trans, alpha, A, x, beta);
% norm(yMAT - yOut)
% Multiply matrix A by row vector x^T with yOut = A^T*x + beta*y
% Trans = 'T';
% alpha = 1;
% A = [1, 2, 3; 4, 5, 6; 7, 8, 9; 10, 11, 12];
% x = double([-1; -2; -3; -4]);
% beta = -2;
% y = [1; -1; 1];
% yMAT = alpha*A'*x + beta*y;
% yOut = dgemv(Trans, alpha, A, x, beta, y);
% norm(yMAT - yOut)
Steps 4 and 5 are now next up: check your code for ways to optimize, keep playing around with it, and maybe even repurpose it into a program of your own for an application.
Now for some practice:
Use what is given here to start and then write your own proper program using function files.
We will be prototyping a function to calculate the outer product of two vectors.
Let 
. Then, the outer product between x and y is the matrix given by
. Then, the outer product between x and y is the matrix given by
Let's code up the algorithm for it using random matrices in a sandbox:
m = 3;
n = 4;
x = randi(10, m, 1);
y = randi(10, n, 1);
AMAT = x*y';
x
y
AMAT
Let's do it using column partition (BLAS 2):
m = 3;
n = 4;
x = randi(10, m, 1);
y = randi(10, n, 1);
AMAT = x*y';
A = zeros(m,n);
for k = 1:n
A(:,k) = y(k,1)*x;
end
err = A - AMAT
By the entries in each column (BLAS 1):
m = 3;
n = 4;
x = randi(10, m, 1);
y = randi(10, n, 1);
AMAT = x*y';
A = zeros(m,n);
for k = 1:n
for i = 1:m
A(i,k) = y(k,1)*x(i,1);
end
end
err = A - AMAT
Through row partition (BLAS 2):
m = 3;
n = 4;
x = randi(10, m, 1);
y = randi(10, n, 1);
AMAT = x*y';
A = zeros(m,n);
for i = 1:m
A(i,:) = x(i,1)*y';
end
err = A - AMAT
By the entries in each row (BLAS 1):
m = 3;
n = 4;
x = randi(10, m, 1);
y = randi(10, n, 1);
AMAT = x*y';
A = zeros(m,n);
for i = 1:m
for k = 1:n
A(i,k) = x(i,1)*y(k,1);
end
end
err = A - AMAT
You can use this testing to write your own function file and then add it to your library.
Let's move on to a rank-one update.
Let 
. Then, a rank-one update of A is given by
. Then, a rank-one update of A is given by
Let's code up the algorithm for it using random matrices in a sandbox.
Here it is by column partition (BLAS 2):
m = 3;
n = 3;
A = rand(m,n);
x = rand(m,1);
y = rand(n,1);
BMAT = A + x*y';
B = zeros(m,n);
for k = 1:n
B(:,k) = A(:,k) + y(k,1)*x;
end
err = B - BMAT
We can make this more memory efficient:
m = 3;
n = 3;
A = rand(m,n);
x = rand(m,1);
y = rand(n,1);
BMAT = A + x*y';
for k = 1:n
A(:,k) = A(:,k) + y(k,1)*x;
end
err = A - BMAT
By the entries in each column (scalar information as opposed to vector information thus BLAS1):
m = 3;
n = 3;
A = rand(m,n);
x = rand(m,1);
y = rand(n,1);
BMAT = A + x*y';
for k = 1:n
for i = 1:m
A(i,k) = A(i,k) + y(k,1)*x(i,1);
end
end
err = A - BMAT
If you don't want to turn either of these into an entire library function, just try writing some specifications.
Here's a start:
where:
α is a scalar
x is vector
y is another vector
A is a matrix
Our function should:
- Take in a scalar α stored in double precision
- Take in an matrix A stored in double precision
- Take input vector x of dimensions
- Take input vector y of dimensions
An error message is produced if...
Good luck with your programming! Have fun!
10. Debugging
When you code, running into bugs is common for even the most senior software engineers. Let's go over some methods you can use to debug your code in MATLAB.
Run and analyze the outputs for certain sections of your code.
- For basic prototyping, you can regularly print outputs to see how things are moving.
- Use robust error checking for better feedback by using error.
- You can divide your code up into sections using Section Break.
- You can run from a specific section to the end of the program using Run to End.
- Use breakpoints to step in and step out of your code.
Let's use the following program as an example. I've moved it to a file named debug1.m:
% Average calculator
data = [5, 10, 15, 20];
n = length(data);
sum = 0;
for i = 1:n+1
sum = sum + data(i);
end
average = sum / n;
fprintf('Average is %.2f\n', average)
As you can see, MATLAB returns an error:
It's not uncommon for bugs to arise from control flow statements. Place a standard breakpoint on the for loop by clicking the line number in the gutter:
Now, let's press Run:
Type dbstep into the Command Window or press Step to move down a line. We can see where we are in Workspace:
We are at i = 1.
Enter dbstep to move to line 8 and then enter it in again to move to the next index:
Keep repeating this process until you find the error:
Let's step once more:
It looks like we found the error. From i = 1 to i = 4, we had no errors. On i = 5, the error appeared. This is because there is no 5th element. Our loop is trying to access an index that does not exist.
Let's fix that. You generally want to avoid editing a file while the program is paused because your changes will not appear unless you rerun the code. We will exit debugging mode and then edit the code:
% Average calculator
data = [5, 10, 15, 20];
n = length(data);
sum = 0;
for i = 1:n % This is where the bug was
sum = sum + data(i);
end
average = sum / n;
fprintf('Average is %.2f\n', average);
Our program is working properly now.
We can also use conditional breakpoints, which allow the debugging process to stop at a certain line under a specific condition. The condition must return a logical scalar value.
This is done by right clicking on the line number where you want the condition evaluated and selecting Set Conditional Breakpoint...:
You then enter your condition here. Here are some example conditions:
- i > length(myArray) is for testing if your index exceeds the array
- denominator == 0 is for testing if you are dividing by 0
- isnan(output(i)) || isinf(output(i)) is for finding invalid or corrupt data
Here is a summary of some of the most important debugging actions:
Action | Function | Description |
Continue  | dbcont | Resume running until the end of the file or another breakpoint |
Step  | dbstep | Execute the current line of code |
Step In  | dbstep in | Execute the current line of code and enter a specific function if called |
Step Out  | dbstep out | After stepping in, run the rest of the called function, and then leave the called function |
Stop  | dbquit | End the debugging session |
Set breakpoint | dbstop | Set a breakpoint at the current line |
Clear breakpoint | dbclear | Clear the breakpoint at the current line |
Here's another helpful function:
help MException
Let's go over a debugging example with the Collatz conjecture. It is one of the most famous unsolved problems in mathematics.
It essentially says that you can take any positive integer and...
- If even, divide it by 2 if even
- If odd, multiply it by 3 and add 1
Repeat with your new number and once you reach a power of 2, you will eventually reach 1.
Can you prove that this will always converge to 1?
Here's an example:
You notice how the calculated integers can go up and down? These sequences are known as hailstone sequences. That's what makes this so difficult to prove true: the unpredictability of numbers. Simplicity can lead to chaos:
max_n = 5000;
iterations = zeros(1, max_n);
% Compute Collatz iteration counts
for n = 1:max_n
count = 0;
current = n;
while current ~= 1
if mod(current, 2) == 0
current = current / 2;
else
current = 3*current + 1;
end
count = count + 1;
end
iterations(n) = count;
end
% Plot
hFig = figure;
plot(1:max_n, iterations, '.', 'MarkerSize', 5);
xlabel('Starting Integer (n)', 'FontSize',14,'FontWeight','bold');
ylabel('Iterations to Reach 1', 'FontSize',14,'FontWeight','bold');
title (['Collatz Iterations for n = 1:' num2str(max_n)], 'FontSize',16,'FontWeight','bold');
grid on;
set(gca, 'FontSize',12, 'LineWidth',1.2);
box on;
How might it be proved false? Let's say you start with 1: the 4, 2, 1 sequence will always repeat. Someone might find a cycle like this that repeats forever.
Paul Erdős, one of the most prolific mathematicians of his time famously said that "mathematics is not yet ready for such problems."
Jeff Lagarias, the mathematician who has perhaps contributed the most to the conjecture in modern day, has warded against mathematicians pursuing this problem.
Terry Tao, the most famous mathematician alive, "proved" that "almost all" starting values will eventually reach 1 in 2019.
Let's move on to the code. This program plots the sequence length versus starting integer. For the first function file, we have:
function sequence = collatz(n)
sequence = n;
next_value = n;
while next_value > 1
if rem(next_value,2)==0
next_value = next_value/2;
else
next_value = 3*next_value+1;
end
sequence = [sequence, next_value];
end
end
For the second we have:
% Plot
function collatzplot(n)
clf
set(gcf,'DoubleBuffer','on')
set(gca,'XScale','linear')
for m = 1:n
plot_seq = collatz(i);
seq_length(m) = length(plot_seq);
line(m, seq_length, 'Marker', '.', 'MarkerSize', 9, 'Color', 'blue');
drawnow
end
end
Running collatzplot(3) returns this plot:
Our algorithm looks sound, so there must be a problem with something in the body of the for loop being used to actually draw the plot. Let's place a breakpoint on this line:
line(m, seq_length, 'Marker', '.', 'MarkerSize', 9, 'Color', 'blue');
For the first index, it looks good:
It looks like the problem starts here. Why does plot_seq have two elements?
And it continues here:
It looks like the entire Collatz sequence is being plotted instead of what we want.
Now that we have isolated the problem, let's look at the code we have and see what we can fix.
It is in the second parameter.
line(m, seq_length, 'Marker', '.', 'MarkerSize', 9, 'Color', 'blue');
It is now fixed.
line(m, seq_length(m), 'Marker', '.', 'MarkerSize', 9, 'Color', 'blue');
Let's plot again:
Now that our function is working, let's try it with a larger integer just for fun.
Lastly, let's learn how to use try/catch statments to silence the effects of errors and potentially give useful debugging info. Take a look at the pseudocode:
try
erroneousOperation
catch
handleError
end
We have a vector with only 3 elements. Let's try to access the 5th:
Let's use them to suppress the runtime errors and instead add our own debugging information:
Further Reading
MATLAB Introductory Guide uses and was inspired by the following sources:
- Jeff Anderson
- Xenophontos — Beginner's Guide to MATLAB
- Houeque — Introduction to MATLAB for Engineering Students
- Pearson — Logical Indexing, NaN, Basic Plotting
- Northwestern ECE — Debugging Example
Good Practices and Style Guides:
- Stepp — CS106X Style Guide
- Johnson — MATLAB Programming Style Guidelines
- Lee — Ten simple rules for documenting scientific software
- Wilson et al — Best Practices for Scientific Computing
- Schlömer — Guidelines for writing clean and fast code in MATLAB