Contents

1. A Few Elementary Calculations

2. Finding Help in Matlab

There are several ways to get help in Matlab:

• help functionname provides direct on-line help on the function whose eaxct name you know. For example help help provides help on the function help.
• lookfor keyword provides a list functions with the string keyword in their description. For example,
  >> lookfor prime  
  isprime                        - True for prime numbers.  
  factor                         - Prime factors.  
  primes                         - Generate list of prime numbers.

  By clicking on isprime, factor, or primes additional help on those functions will be listed.

• helpwin provides the click and navigate help. A new help windows opens. One can also select help choices from the Help menu.
• helpdesk is the web browser-based help, including printable (PDF) documentation on the Web.
• demo provides the demonstration programs on Matlab.

Additionally, one can find detailed online help on Matlab at:

http://www.mathworks.com/access/helpdesk/help/techdoc/index.html

3. Working with Arrays of Numbers

4. Creating Script Files and Simple Plots

A script file is an external file that contains a sequence of Matlab statements. By typing the filename, subsequent Matlab input is obtained from the file. Script files have a filename extension of .m and are often called M-files.

Select New from File menu (alternatively, use any text editor) and type the following lines into this file:

Lines starting with % sign are ignored; they are interpreted as comment lines by Matlab.

Select Save As… from File menu and type the name of the document as, e.g., ellipse_parabola.m. Alternatively, if you are using an external editor, use a suitable command to save the document in the file with .m extension.

Once the file is placed in the current working directory of Matlab, type the following command in the command window to execute ellipse_parabola.m script file.

The above script file also illustrates how to generate plot $x$ vs $y$ from the generated $x$ and $y$ coordinates of 100 equidistant points between $0$ and $2\pi$. For example, in the above script change the entry plot(x,y,t,z) to plot(x,y,’--’,t,z,’o’) and see the difference in graphs drawings. For more information, check Matlab help on plot command and its options by invoking help plot in the command window of Matlab. For print options enter help print in the command window of Matlab.

5. Creating and Executing a Function File

Function file is a script file (M-file) that adds a function definition to Matlab’s list of functions. Edit ellipse_parabola.m script file (listed in Section 4 to look like the following:

Save it in circle_f1.m file (remember to place it in the working directory of Matlab) and try it out.

And here is slightly different version of circle_f1.m function file, called circle_f2.m This function file, in addition to drawing a circle with radius $r$, computes the area inside this circle together with the outputs for x and y coordinates.

Save it in circle_f2.m file (remember to place it in the working directory of Matlab) and execute it in one of the following ways:

Finally, suppose you need to evaluate the following function many times at different values of $x$ during your Matlab session: $f\left(x\right)=exp\left(-0.3x\right)sin\left(30x\right)+\frac{x^{10}}{1+x^{10}}$. You can create an inline function and pass its name to other functions:

In the case you want to use function $f\left(x\right)=exp\left(-0.3x\right)sin\left(30x\right)+\frac{x^{10}}{1+x^{10}}$ in future Matlab sessions, create the following function file:

and save it in my_f.m file (remember to place it in the working directory of Matlab). You can use it in many ways:

And here is Matlab function that for a given $x\in ℝ$ and given $n\in \mathbb{N}$ computes

Save the script in expseries.m file (remember to place it in the working directory of Matlab) and try it out. Alternatively, you can also download it from:

http://www.csun.edu/~hcmth008/matlab_tutorial/expseries.m

6. Saving/Loading Data, Recording a Session, Clearing Variables, and More on Plotting

In addition to plot command, introduced in Section 4, Matlab has many plotting utilities. Here are sample uses of two of them.

For other plotting utilities check Matlab help on the following commands:
ezcontour, ezcontourf, ezmesh, ezmeshc, ezplot3, ezpolar, ezsurf, ezsurfc.

7. Matrices and Vectors

The command linspace(a,b,n) creates a linearly spaced (row) vector of length $n$ from $a$ to $b$. Note that u=linspace(a,b,n) is the same as the command u=a:(b-a)/(n-1):b.

Also, logspace(a,b,n) generates a logarithmically spaced (row) vector of length $n$ from $10^a$ to $10^b$. For example,

Finally, logspace(a,b,n) is the same as 10.^(linspace(a,b,n) where .^ is term by term operator of exponentiation (see Section 7.3).

There are many ways to create special vectors, for example, vectors of zeros or ones of specific length:

7.1. Appending or Deleting a Row or Column.

A row/column can be easily appended or deleted, provided the row/column has the same length as the length of the rows/columns of the existing matrix.

7.2. Matrix Operations.

In Matlab, as in Pascal, Fortran, or C, matrix product is just C=A*B, where $A$ is an $m\times n$ matrix and $B$ is $n\times k$ matrix, and the resulting matrix $C$ has the size $m\times k$.

• A+B and A-B are valid if the matrices $A$ and $B$ have the same size (i.e., the same number of rows and the same number of columns).
• A*B is valid if $A$’s number of columns equals $B$’s number of rows. The command A^2, which equals A*A, makes sense if $A$ is a square matrix.
• X=A/B (right division) is valid if $A$ and $B$ have the same number of columns and it denotes the solution to the matrix equation $XA=B$. In particular, if $A$ and $B$ are square matrices, $A∕B\equiv A{B}^{-1}$, if $B^{-1}$ exists.
• X=A\B (left division) is valid if $A$ and $B$ have the same number of rows and it denotes the solution to the matrix equation $AX=B$. In particular, if $A$ and $B$ are square matrices, $A\setminus B\equiv A^{-1}B$, if $A^{-1}$ exists.

If $A$ is $m\times p$ and $B$ is $p\times n$, their product $C=AB$ is $m\times n$. Matlab defines their product by C = A*B, though this product can be also defined using Matlab for loops (see Section 9) and colon notation considered in the Section 7.1. However, the internally defined vectorized form of the product A*B is more efficient; in general, such vectorizations are strongly recommended, whenever possible.

The matrix A=magic(3) belongs to the family of so called magic squares. M = magic(n) returns an $n\times n$ matrix ($n\ge 3$) constructed from the integers 1 through $n^2$ such that the $n$ numbers in all rows, all columns, and both diagonals sum to the same constant. The constant sum in every row, column and diagonal is called the magic constant or magic sum $M$. The magic constants of a magic square depend only on $n$ and have the values

(See, for example, http://en.wikipedia.org/wiki/Magic_square .)

The commands below check that the sums in all rows, columns, and in main diagonals of A=magic(5) are indeed the same and equal to 65.

For information on pascal(n) matrices see Section 12.1.

The scalar (dot) product of two vectors or multiplication of a matrix with a vector can be done easily in Matlab.

7.3. Array Operations.

Array operations are done on an element-by-element basis for matrices/vectors of the same size:

• .* denotes element-by-element multiplication;
• ./ denotes element-by-element left division;
• .\ denotes element-by-element right division;
• .^ denotes element-by-element exponentiation.

7.4. Matrix Functions.

• expm(A) computes the exponential of a square matrix $A$, i.e., $e^A$. This operation should be distinguished from exp(A) which computes the element-by-element exponential of $A$;
• logm(A) computes $log\left(A\right)$ of square matrix $A$. $log\left(A\right)$ is the matrix such that $A=e^{log\left(A\right)}$;
• sqrtm(A) computes $\sqrt{A}$ of square matrix. $\sqrt{A}$ is the matrix such that $A=\left(\sqrt{A}\right)\ast \left(\sqrt{A}\right)$.

For more information on these functions see also Section 12.8.6.

8. Character Strings

S = ’Any Characters’ creates a character array, or string. All character strings are entered using single right-quote characters. The string is actually a vector whose components are the numeric codes for the characters (the first 127 codes are ASCII). The length of S is the number of characters. A quotation within the string is indicated by two single right-quotes.

S = [S1 S2 ...] concatenates character arrays S1, S2, etc. into a new character array, S.

S = strcat(S1, S2, ...) concatenates S1, S2, etc., which can be character arrays or Cell Arrays of Strings.

Trailing spaces in strcat character array inputs are ignored and do not appear in the output. This is not true for strcat inputs that are cell arrays of strings. Use the S = [S1 S2 ...] concatenation syntax, to preserve trailing spaces.

A column vector of text objects is created with C=[S1; S2; ...], where each text string must have the same number of characters.

The command char(S1,S2,...) converts strings to a matrix. It puts each string argument S1, S1, etc., in a row and creates a string matrix by padding each row with the appropriate number of blanks; in other words, the command char does not require the strings S1, S2, …to have the same number of characters. One can also use strvcat and str2mat to create strings arrays from strings arguments, however, char does not ignore empty strings but strvcat does.

8.1. eval Function.

eval executes string containing Matlab expression. For example, the command

evaluates $y=exp\left(3\right)cos\left(\pi ∕4\right)$ and is equivalent to typing $y=exp\left(3\right)cos\left(\pi ∕4\right)$ on the command prompt. eval function is often used to load or create sequentially numbered data files. Here is a hypothetical example:

These are the strings being evaluated:

9. for and while Loops and if Conditionals

for x=initval:endval, statements end command repeatedly executes one or more Matlab statements in a loop. Loop counter variable x is initialized to value initval at the start of the first pass through the loop, and automatically increments by 1 each time through the loop. The program makes repeated passes through statements until either x has incremented to the value endval, or Matlab encounters a break, or return instruction,

while expression, statements, end command repeatedly executes one or more Matlab statements in a loop, continuing until expression no longer holds true or until Matlab encounters a break, or return instruction.

if conditionally executes statements. The general form of the if command is

The else and elseif parts are optional.

Example 1.

The graph of function $f\left(x\right)=3x^2-e^x$

indicates that $f\left(x\right)=3x^2-e^x$ has two positive zeros.

The computation of these two positive roots can be achieved by the so called fixed-point iteration scheme. For example, the smaller positive root $p$ is the fixed point of $g_1\left(x\right)=\sqrt{\left(\frac{1}{3}\cdot e^x\right)}$. defined on $\left[0,1\right]$, i.e., such that

The larger positive root is the fixed point $c$ of $g_2\left(x\right)=ln\left(3x^2\right)$ defined on $\left[3,4\right]$, i.e., such that

These fixed points (and ultimately the positive roots of $f\left(x\right)=3x^2-e^x$) can be computed accurately by running the so called fixed-point iteration scheme in the form:

And here is the actual Matlab session for the above fixed-point iteration scheme based on $g_1$:

and here fixed-point iteration scheme based on $g_2$:

The accurate value of the root of $f\left(x\right)=3x^2-e^x$ in $\left[0,1\right]$ (to within $10^{-10}$) is $p=0.9100075725$.

The accurate value of the root of $f\left(x\right)=3x^2-e^x$ in $\left[3,4\right]$ (to within $10^{-10}$) is $p=3.733079029$.

Note: It is interesting to know that the fixed-point iteration based on $g_1$ does not converge to the larger positive root of $f\left(x\right)=3x^2-e^x$. Similarly, the fixed-point iteration based on $g_2$ does not converge to the smaller positive root of $f\left(x\right)=3x^2-e^x$.

Example 2.

The Euler numerical method for approximating the solution to the differential equation,

has the form

In the above scheme, $y_n$, for $n=1,2,\dots$, approximates true solution at $x_n$, i.e., $y\left(x_n\right)$.

The initial value problem

has the exact solution

We will use the Euler method to approximate this solution. We will also compute the differences between the exact and approximate solutions.

>> error = abs(Y_exact-y)’   error =            0       0.0587       0.0793       0.0794       0.0696       0.0560       0.0418       0.0289       0.0180       0.0092       0.0025       0.0024       0.0058       0.0080       0.0093       0.0099       0.0100       0.0097       0.0092       0.0086       0.0079 Figure 3: Plots of the exact solution (solid line) and the approximate solution (dashed line) on the interval $\left[0,2\right]$.

10. Interactive Input

The commands input, keyboard, menu, and pause provide interactive user input and can be used inside a script or function file.

keyboard command, when placed in a script or function file, stops execution of the file and gives control to the user’s keyboard. The special status is indicated by a K appearing before the prompt. Variables may be examined or changed - all Matlab commands are valid. The keyboard mode is terminated by executing the command return and pressing the return key. Control returns to the invoking script/function file.

Command menu generates a menu of choices for user input. On most graphics terminals menu will display the menu-items as push buttons in a figure window, otherwise they will be given as a numbered list in the command window. For example,

creates a figure with buttons labeled Red, Blue and Green. The button clicked by the user is returned as choice (i.e. choice = 2 implies that the user selected Blue).

After input is accepted, the dialog box closes, returning the output in choice. You can use choice to control the color of a graph:

pause commands waits for user response. pause(n) pauses for n seconds before continuing.

pause causes a procedure to stop and wait for the user to strike any key before continuing.

pause off indicates that any subsequent pause or pause(n) commands should not actually pause. This allows normally interactive scripts to run unattended.

pause on (default setting) indicates that subsequent pause commands should pause.

pause query returns the current state of pause, either on or off.

11. Input/Output Commands

Matlab provides many C-language file I/O functions for reading and writing binary and text files. The following five I/O functions: fopen, fclose, fgetl, fscanf, and fprintf are frequently used.

• fid = fopen(filename) opens or creates the filename. The filename argument is a string enclosed in single quotation marks.
  Output value fid is a scalar Matlab integer that you use as a file identifier for all subsequent low-level file input/output routines.
• fclose(fid) closes the specified file if it is open, returning 0 if successful and -1 if unsuccessful. fid is an integer file identifier associated with an open file.
  fclose(’all’) closes all open files.
• tline = fgetl(fid) returns the next line of the file associated with the file identifier fid. The returned string tline does not include the newline characters with the text line. To obtain the newline characters, use fgets function.
• A = fscanf(fid, format) reads data from the text file specified by fid, converts it according to the specified format string, and returns it in matrix A. Argument fid is an integer file identifier obtained from fopen. format is a string specifying the format of the data to be read.
• count = fprintf(fid, format, A, ...) formats the data in the real part of matrix A (and in any additional matrix arguments) under control of the specified format string, and writes it to the file associated with file identifier fid. count returns a count of the number of bytes written.
  Argument fid is an integer file identifier obtained from fopen. It can also be 1 for standard output (the screen) or 2 for standard error. Omitting fid causes output to appear on the screen. For example, fprintf(’A unit circle has circumference %g radians.\n’, 2*pi) displays a line on the screen:
  A unit circle has circumference 6.283186 radians.

The example below reads and displays the first three lines from the file circle_f2.m (see, Section 5) and available from: http://www.csun.edu/~hcmth008/matlab_tutorial/circle_f2.m
It also tests for end-of-file and then closes circle_f2.m file.

The value of variable is_end_of_file = 0 indicates that the end of circle_f2.m file was not reached.
Changing name of the file and replacing i<=3 by i<=n will make the above script read and display the first n lines. The entire file will be displayed if the number of lines in the file is less or equal to n.

The next script reads and displays every line from the file my_f.m (created in Section 5) and available from: http://www.csun.edu/~hcmth008/matlab_tutorial/my_f.m
It also tests for end-of-file and then closes the file my_f.m. No a priori information about the number of lines in the file is needed in the script.

The value of variable is_end_of_file = 1 indicates that the end of my_f.m file was reached.

The next example is more advanced. It shows how to read and display lines from given file in a non-sequential order. It uses eval and int2str commands for a customized display. The script below reads and displays the 4th, 5th, 8th, and 9th lines from circle_f2.m file. It also appends the numbers to the above lines. No a priori information about the number of lines in the file is needed in the script.

The next examples show how to use fprintf and fscanf commands. First, create a text file called exp.txt containing a short table of the exponential function. (On Windows platforms, it is recommended that you use fopen with the mode set to ’wt’ to open a text file for writing.)

Now, examine the contents of exp.txt with type filename command.

One can now read this file back into a two-column Matlab matrix with the script:

Comparing it with exp.txt and long format outputs we observe no loss of accuracy.

As the final example, start with a file temp.dat that contains temperature readings

File temp.dat is available from: http://www.csun.edu/~hcmth008/matlab_tutorial/temp.dat

Open the file using fopen and read it with fscanf. If you include ordinary characters (such as the degrees character ${}^{\circ }$ (ASCII 176) and Fahrenheit (F) symbols used here) in the conversion string, fscanf skips over those characters when reading the string. This allows for further processing (see below). The display on the right is given for a comparison.

12. Linear Algebra

12.1. Solving Linear Square Systems.

Write the system of linear algebraic equations

in matrix form

with

Please note that the number of columns of $A$ equals to the number of rows of $X$, and thus the matrix multiplication of $A$ with $X$ is defined.

Enter the matrix $A$ and vector $b$, and solve for vector $X$ with A\b.

The matrix $A$ belongs to the family of so called Pascal matrices and can be entered by the command A = pascal(4).

In general, A = pascal(n) returns the Pascal matrix of order n: a symmetric positive definite matrix with integer entries taken from Pascal’s triangle that determines the coefficients which arise in binomial expansions. See, for example, http://en.wikipedia.org/wiki/Pascal’s_triangle .

12.1.1. Singular Coefficient Matrix. A square matrix $A$ is singular if it does not have linearly independent columns. $A$ is singular if and only if the determinant of $A$ is zero (see Section 12.3). If $A$ is singular, the equation $AX=b$ either has no solution or it has infinite many solutions. Actually, one can show that $AX=b$ system has a solution if and only if $b\perp N$ (i.e, the scalar product of $b$ and $N$ is zero) for all $N$ such that $A^{\prime }N=0$, where $A^{\prime }$ is the transpose (the complex conjugate transpose, if $A$ has complex entries) of $A$, i.e., the matrix where the original rows and columns are interchanged. The matrix

is singular, as easy verified by typing

For b = [5;2;12], the equation $AX=b$ has a solution given by

For more information on pinv(A), see Section 12.4. In this case, command A\b returns an error condition.

One can also verify that b = [5;2;12] is perpendicular to any solution of the homogeneous system $A^{\prime }N=0$:

and

This verifies that $b\perp N$.

The general (complete) solution to $AY=b$ is given by

where $Z$ solves the homogeneous system $AZ=0$ and $c$ is any constant (scalar), as verified by

On the other hand, for b = [1;2;3], $b\perp ̸N$. (indeed, dot(N,b)=-0.4082) and the system $AX=b$ does not have an exact solution.

In this case, pinv(A)*b returns a least squares solution $X$, i.e.,

where for $v=\left(v_1,v_2,\dots ,v_n\right)$, $\parallel v\parallel ={\left(\sum _{i=1}^n|v_i{|}^2\right)}^{1∕2}$ is the so called norm of vector $v$. Matlab command for the norm of vector v is norm(v). We have

12.2. Gaussian Elimination.

One of the techniques learned in linear algebra courses for solving a system of linear equations is Gaussian elimination. One forms the augmented matrix that contains both the coefficient matrix $A$ and vector $b$, and then, using Gauss-Jordan reduction procedure, the augmented matrix is transformed into the so called row reduced echelon form. Matlab has built-in function that does precisely this reduction. For A = pascal(4) and b = [1; -1; 2; 3], considered in Section 12.1, the commands C=[A b] and rref(C) calculate the augmented matrix and the row reduced echelon matrix, respectively.

Since matrix A=pascal(4) is not singular, the last column of rref(C) (One could also use rref([A b]) command here.) provides the solution X already obtained in Section 12.1 by using X = A\b command.

Here is another use of rref command. One can determine whether $AX=b$ has an exact solution by finding the row reduced echelon form of the augmented matrix [A b]. For example, enter

Since the bottom row contains all zeros except for the last entry, the equation $AX=b$ does not have a solution. Indeed, the above row reduced echelon form is equivalent to the following system of equations:

As indicated in Section 12.4, in this case, pinv(A)*b command returns a least squares solution.

12.3. Inverses, Determinants, and Cramer’s Rule.

For a square and nonsingular matrix $A$, the equations $AX=I$ and $XA=I$ (here $I$ is an identity matrix), have the same solution, $X$. This solution is called the inverse of $A$ and is denoted by $A^{-1}$. Matlab function inv(A) computes the inverse of $A$.

The inverse of A=pascal(4) is

and thus, the solution of the system $AX=b$ with b=[1;0;2;0] is

The determinant is a special number associated with any square matrix $A$ and is denoted by $\left|A\right|$ or $det\left(A\right)$. The meaning of a determinant is a scale factor for measure when the matrix is regarded as a linear transformation. Thus a $2\times 2$ matrix, with determinant 3, when applied to a set of points with finite area will transform those points into a set with 3 times the area.

The matrix $A=\left[\begin{array}{cc}\hfill a\hfill & \hfill b\hfill \\ \hfill c\hfill & \hfill d\hfill \end{array}\right]$ has the determinant $det\left(A\right)=ad-bc$ and the geometric interpretation presented in Figure 4. Figure 4: The area of the parallelogram is the absolute value of the determinant of the matrix formed by the vectors representing the parallelogram’s sides.

The determinant of

and the Rule of Sarrus is a mnemonic for this formula:

Figure 5: The determinant of a $3\times 3$ matrix can be calculated by its diagonals. Figure 6: The volume of this parallelepiped is the absolute value of the determinant of the matrix formed by the rows r1, r2, and r3.

The determinant of $n\times n$ matrix $A={\left(a_{ij}\right)}_{i,j=1,2,\dots n}$ ($n\ge 1$) can be defined by the Leibniz formula:

Here the sum is computed over all permutations $\sigma$ of the numbers $\left\{1,2,\dots ,n\right\}$. A permutation is a function that reorders this set of integers. The set of such permutations is denoted by $S_n$. The sign function of permutations in the permutation group $S_n$ returns $+1$ and $-1$ for even and odd permutations, respectively. For example, for $n=4$ and $\sigma =\left(1,4,3,2\right)$, $\text{sign}\left(\sigma \right)=-1$ (one pair switch) and $\prod _{i=1}^4{a}_{i\sigma \left(i\right)}=a_{11}{a}_{24}{a}_{33}{a}_{42}$.

Furthermore, for $A=\left[\begin{array}{cc}\hfill a_{11}\hfill & \hfill a_{12}\hfill \\ \hfill a_{21}\hfill & \hfill a_{22}\hfill \end{array}\right]$, $S_2=\left\{{\sigma }_1,{\sigma }_2\right\}$, with ${\sigma }_1=\left(1,2\right)$, ${\sigma }_2=\left(2,1\right)$, $\text{sign}\left({\sigma }_1\right)=+1$, $\text{sign}\left({\sigma }_2\right)=-1$, and we have

Matlab function det(A) computes the determinant of a square matrix A. For A = pascal(4), we have

Actually, det(pascal(n)) is equal to 1 for any $n\ge 1$.

A square matrix $A$ has the inverse $A^{-1}$ if and only if $det\left(A\right)\ne 0$. In this case, $det\left(A^{-1}\right)=1∕det\left(A\right)$.

12.3.1. Cramer’s rule. For square matrices $A$ that are invertible, the system of equations $AX=b$ with a given $n$-dimensional column vector $b$ can be solved by using so called Cramer’s rule. This rule states that the solution $X={\left(x_1,x_2,\dots ,x_n\right)}^T$ is obtained by the formula:

where $A_i$ is the matrix formed by replacing the ith column of $A$ by the column vector $b$.

Here is Matlab code for Cramer’s rule applied to A=magic(5) and b=[1;2;3;4;5]

Check that the backslash operator A\b and the command inv(A)*b provide the same solution.

Note: For large size matrices, the use of Cramer’s rule and of the inverse is not computationally efficient as compared to the backslash operator.

12.4. Pseudoinverses or Generalized Inverses.

The pseudoinverse $A^{+}$ f an $m\times n$ matrix $A$ is a generalization of the inverse matrix. The terms Moore-Penrose pseudoinverse or generalized inverse are also used as synonyms for pseudoinverse. Pseudoinverse is a unique $n\times m$ matrix $A^{+}$ satisfying all of the following four conditions:

1. $A{A}^{+}A=A$;
2. $A^{+}A{A}^{+}=A^{+}$;
3. ${\left(A{A}^{+}\right)}^{\prime }=A{A}^{+}$ (i.e.,$A{A}^{+}$ is Hermitian);
4. ${\left(A^{+}A\right)}^{\prime }=A^{+}A$ (i.e.,$A^{+}A$ is Hermitian).

Here $A^{\prime }$ is the Hermitian transpose (also called conjugate transpose) of a matrix $A$. For matrices whose elements are real numbers instead of complex numbers, $A^{\prime }=A^T$, In Matlab notation, Hermitian transpose of A is A’.

Pseudoinverses satisfy many properties. Here are some of them:

• The pseudoinverse exists and is unique for any matrix $A$.
• If $A$ has real entries, then so does $A^{+}$; if $A$ has rational entries, then so does $A^{+}$.
• If the matrix $A$ is invertible, the pseudoinverse and the inverse coincide, i.e., $A^{+}=A^{-1}$.
• The pseudoinverse of a zero matrix is its transpose.
• The pseudoinverse of the pseudoinverse is the original matrix, i.e., ${\left(A^{+}\right)}^{+}=A$.
• Pseudoinversion commutes with transposition, conjugation, and taking the conjugate transpose.
• The pseudoinverse of a scalar multiple of $A$ is the reciprocal multiple of $A^{+}$, i.e., ${\left(\alpha A\right)}^{+}={\alpha }^{-1}{A}^{+}$ for $\alpha \ne 0$.

For more information go to the webpage: http://en.wikipedia.org/wiki/Moore-Penrose_pseudoinverse

A common use of the pseudoinverse is to compute a best fit (least squares) solution to a system of linear equations that lacks a unique solution. In general, a given linear system $AX=b$ may not have a solution. Even if there exists such a solution, it may not be unique. One can however ask for a vector $X$ that brings $AX$ as close as possible to vector $b$, i.e., a vector that minimizes the Euclidean norm

where for $v=\left(v_1,v_2,\dots ,v_n\right)$, $\parallel v\parallel ={\left(\sum _{i=1}^n|v_i{|}^2\right)}^{1∕2}$. Matlab command for the norm of vector v is norm(v).

If there are several such vectors $X$, one could ask for the one among them with the smallest Euclidean norm. Thus formulated, the problem has a unique solution, given by the pseudoinverse:

Matlab function that computes the pseudoinverse of a matrix A is pinv(A). If for given A and $$b, the least square solution X is unique, it can be computed by entering either X=pinv(A)*b or X=A\b. If there exist many least square solutions (for the so called rank-deficient matrices), X=pinv(A)*b computes vector X with the smallest norm norm(X), while X=A\b computes a vector X that has at most $r$ nonzero components. Here, $r$ is the rank of matrix A equal to rank(A).

Sections 12.5.1 and 12.5.2 provide many examples for the use of pinv(A).

12.5. Underdetermined and Overdetermined Linear Systems.

12.5.1. Underdetermined linear systems. Underdetermined linear systems involve more unknowns than equations. Such systems are studied in linear programming, especially when accompanied by additional constraints. A solution, when exists, is never unique. For A=[0 1 1 0; 1 1 0 1; 1 2 1 1] and b=[1;2;1], the corresponding system has three equations and four unknowns:

The row reduced echelon form shows that the system has no exact solutions. In this case, least square solutions $X$ such that

can be obtained by entering A\b or pinv(A)*b (see Section 12.4):

We observe that there are many vectors Z (here, X and Y) that minimize norm(A*Z-b). This is due to the fact that matrix $A$ has only 2 linearly independent columns (so called rank-deficient system, while the full rank would be 3), the least square problem is not unique.

Consider now the following two equations with four unknowns:

Enter A=[6 8 7 3;3 5 4 1] and b=[1;2], and calculate rref([A b]):

Thus, the system has an exact solution that can be obtained by entering X=A\b or W=pinv(A)*b:

Both X and W are exact solutions:

The general (complete) solution to $AY=b$ is given by

where $c_1$ and $c_2$ are any constants (scalars), and $Z_1$ and $Z_2$ are fundamental solutions of the system $AZ=0$, as verified by

The other solution, W=pinv(A)*b, is also of the form $c_1{Z}_1+c_2{Z}_2+X$, for some constants $c_1$ and $c_2$. Indeed, the constants $c_1$ and $c_2$ are the solutions of the following linear system D*C=W-X with C=[c1;c2], as verified by

12.5.2. Overdetermined linear systems. Overdetermined linear systems involve more equations than unknowns. For example, they are encountered in various kinds of curve fitting to experimental data (see Section 13). Enter the following data into Matlab:

t=[0;0.25;0.5;;0.75;1.00;1.25;1.50;1.75;2.00;2.25;2.50];
y=sin(t)+rand(size(t))/10;

We will model the above randomly corrupted sine function data with a quadratic polynomial function.

In other words, we want to approximate the vector y by a linear combination of three other vectors:

The above system has eleven equations and three unknowns: $c_1$, $c_2$, and $c_3$. It is represented 11-by-3 matrix A and b=y. The row reduced echelon form of A shows the above system does not have an exact solution:

Since rank(A)=3, both the backslash operator (i.e., A\b) and pinv(A)*b provide the same (unique) least square solution:

In other words, the least squares fit to the data is

The following commands evaluate the model at regularly spaced increments in $t$ (finer than the original spacing) and then plot the result including the original data.

For more information on polyfit, see Section 13.1 on curve fitting and interpolation.

Note: Due to randomness of y = sin(t) + rand(size(t))/10, your results will differ slightly from the above interpolating polynomial and the graph.

And here is another attempt at fitting randomly corrupted sine function data into

And here is another attempt at fitting randomly corrupted sine function data into