Episode #125 of the Stack Overflow podcast is here. We talk Tilde Club and mechanical keyboards. Listen now

4 added 485 characters in body

### Big idea: solve linear systems rather than compute matrix inverse and multiply

Computing the inverse of a large matrix is slow and imprecise. If you have a huge matrix, don't compute the inverse. Instead:

• Observe $$A^{-1} \mathbf{b}$$ is the solution to $$A\mathbf{x} = \mathbf{b}$$. More generally, $$A^{-1}B$$ is the solution to $$AX=B$$. Anytime you see, $$A^{-1}B$$ you should think, "Hey, maybe I should solve linear systems rather than compute the matrix inverse?" (Note: for small problems, it probably doesn't matter.)

• Another idea: solving triangular systems is super fast! $$\left[\begin{array}{ccc} r_{11}&r_{12} & r_{13}\\0&r_{22}&r_{23}\\0&0&r_{33} \end{array} \right] \left[ \begin{array}{c}x_1\\x_2\\x_3\end{array} \right] = \left[ \begin{array}{c}b_1\\b_2\\b_3\end{array} \right]$$ Notice you immediately have $$x_{3} = b_3 / r_{33}$$. Then substitute $$x_3$$ and second row $$x_2r_{22}+x_3r_{23}=b_2$$ immediately gives you $$x_2$$. Computers can solve triangular systems at lightning speed applying this backsolve algorithm. (Don't reinvent the wheel, call Linear Algebra libraries to do this.) This is what Gordon is getting at when advocating taking a Cholesky decomposition on $$A$$. Once you have the Cholesky Decomposition $$A=RR'$$ where $$R$$ is lower triangular, you can solve $$A\mathbf{x}=\mathbf{b}$$ with two back substitutions, which is nearly free.

$$\hat{\boldsymbol\beta}=(\mathbf{X^{'}C^{-1}X})^{-1}\mathbf{X^{'}C^{-1}y}$$

If you Cholesky decomposition for $$C$$ such that $$C = RR'$$ and $$C^{-1} = R'^{-1}R^{-1}$$ then:

$$\hat{\boldsymbol\beta}=(X^{'}R'^{-1}R^{-1} X)^{-1}X^{'}R'^{-1}R^{-1}y$$

So you want to solve:

• $$R^{-1}X$$ with a backsolve
• $$R^{-1}y$$ with a backsolve
• and then one final linear system

### Practical coding way forward

Since there's already an answer in R, I'll give some code in Matlab. In Matlab, there's the operator \ which basically examines your matrix and does the right thing to solve a linear system fast.

C_inv_X = C \ X;                     % Solves C^{-1}X
bhat = (X'*C_inv_X) \ (C_inv_X'*y);  % Final linear system solve


You could also:

R       = chol(C);  %Careful... does your library give you an UT or LT?
C_inv_X                    %chol gives you an upper triangular in Matlab
Rp_inv_X = R\(R'\X);R'\X; %mldivide auto-detects triangular system
% in other environments, you may have to declare it to
% get full speed
Rp_inv_y = R'\y;
bhat    = (X'*C_inv_XRp_inv_X'*Rp_inv_X) \ (C_inv_X'*yRp_inv_X'*Rp_inv_y);  % Final linear system solve


And then you have the triangular matrix $$R$$ around to use later.

### Big idea: solve linear systems rather than compute matrix inverse and multiply

Computing the inverse of a large matrix is slow and imprecise. If you have a huge matrix, don't compute the inverse. Instead:

• Observe $$A^{-1} \mathbf{b}$$ is the solution to $$A\mathbf{x} = \mathbf{b}$$. More generally, $$A^{-1}B$$ is the solution to $$AX=B$$. Anytime you see, $$A^{-1}B$$ you should think, "Hey, maybe I should solve linear systems rather than compute the matrix inverse?" (Note: for small problems, it probably doesn't matter.)

• Another idea: solving triangular systems is super fast! $$\left[\begin{array}{ccc} r_{11}&r_{12} & r_{13}\\0&r_{22}&r_{23}\\0&0&r_{33} \end{array} \right] \left[ \begin{array}{c}x_1\\x_2\\x_3\end{array} \right] = \left[ \begin{array}{c}b_1\\b_2\\b_3\end{array} \right]$$ Notice you immediately have $$x_{3} = b_3 / r_{33}$$. Then substitute $$x_3$$ and second row $$x_2r_{22}+x_3r_{23}=b_2$$ immediately gives you $$x_2$$. Computers can solve triangular systems at lightning speed applying this backsolve algorithm. (Don't reinvent the wheel, call Linear Algebra libraries to do this.) This is what Gordon is getting at when advocating taking a Cholesky decomposition on $$A$$. Once you have the Cholesky Decomposition $$A=RR'$$ where $$R$$ is lower triangular, you can solve $$A\mathbf{x}=\mathbf{b}$$ with two back substitutions, which is nearly free.

### Practical coding way forward

Since there's already an answer in R, I'll give some code in Matlab. In Matlab, there's the operator \ which basically examines your matrix and does the right thing to solve a linear system fast.

C_inv_X = C \ X;                     % Solves C^{-1}X
bhat = (X'*C_inv_X) \ (C_inv_X'*y);  % Final linear system solve


You could also:

R       = chol(C);  %Careful... does your library give you an UT or LT?
C_inv_X = R\(R'\X); %mldivide auto-detects triangular system
% in other environments, you may have to declare it to
% get full speed
bhat    = (X'*C_inv_X) \ (C_inv_X'*y);  % Final linear system solve


And then you have the triangular matrix $$R$$ around to use later.

### Big idea: solve linear systems rather than compute matrix inverse and multiply

Computing the inverse of a large matrix is slow and imprecise. If you have a huge matrix, don't compute the inverse. Instead:

• Observe $$A^{-1} \mathbf{b}$$ is the solution to $$A\mathbf{x} = \mathbf{b}$$. More generally, $$A^{-1}B$$ is the solution to $$AX=B$$. Anytime you see, $$A^{-1}B$$ you should think, "Hey, maybe I should solve linear systems rather than compute the matrix inverse?" (Note: for small problems, it probably doesn't matter.)

• Another idea: solving triangular systems is super fast! $$\left[\begin{array}{ccc} r_{11}&r_{12} & r_{13}\\0&r_{22}&r_{23}\\0&0&r_{33} \end{array} \right] \left[ \begin{array}{c}x_1\\x_2\\x_3\end{array} \right] = \left[ \begin{array}{c}b_1\\b_2\\b_3\end{array} \right]$$ Notice you immediately have $$x_{3} = b_3 / r_{33}$$. Then substitute $$x_3$$ and second row $$x_2r_{22}+x_3r_{23}=b_2$$ immediately gives you $$x_2$$. Computers can solve triangular systems at lightning speed applying this backsolve algorithm. (Don't reinvent the wheel, call Linear Algebra libraries to do this.) This is what Gordon is getting at when advocating taking a Cholesky decomposition on $$A$$. Once you have the Cholesky Decomposition $$A=RR'$$ where $$R$$ is lower triangular, you can solve $$A\mathbf{x}=\mathbf{b}$$ with two back substitutions, which is nearly free.

$$\hat{\boldsymbol\beta}=(\mathbf{X^{'}C^{-1}X})^{-1}\mathbf{X^{'}C^{-1}y}$$

If you Cholesky decomposition for $$C$$ such that $$C = RR'$$ and $$C^{-1} = R'^{-1}R^{-1}$$ then:

$$\hat{\boldsymbol\beta}=(X^{'}R'^{-1}R^{-1} X)^{-1}X^{'}R'^{-1}R^{-1}y$$

So you want to solve:

• $$R^{-1}X$$ with a backsolve
• $$R^{-1}y$$ with a backsolve
• and then one final linear system

### Practical coding way forward

Since there's already an answer in R, I'll give some code in Matlab. In Matlab, there's the operator \ which basically examines your matrix and does the right thing to solve a linear system fast.

C_inv_X = C \ X;                     % Solves C^{-1}X
bhat = (X'*C_inv_X) \ (C_inv_X'*y);  % Final linear system solve


You could also:

R       = chol(C);  %Careful... does your library give you an UT or LT?
%chol gives you an upper triangular in Matlab
Rp_inv_X = R'\X; %mldivide auto-detects triangular system
% in other environments, you may have to declare it to
% get full speed
Rp_inv_y = R'\y;
bhat    = (Rp_inv_X'*Rp_inv_X) \ (Rp_inv_X'*Rp_inv_y);  % Final linear system solve


And then you have the triangular matrix $$R$$ around to use later.

3 added 122 characters in body

### Big idea: solve linear systems rather than compute matrix inverse and multiply

Computing the inverse of a large matrix is slow and imprecise. If you have a huge matrix, don't do itcompute the inverse. Instead:

• Observe $$A^{-1} \mathbf{b}$$ is the solution to $$A\mathbf{x} = \mathbf{b}$$. More generally, $$A^{-1}B$$ is the solution to $$AX=B$$. Anytime you see, $$A^{-1}B$$ you should think, "Hey, maybe I should solve linear systems rather than compute the matrix inverse?" (Note: for small problems, it probably doesn't matter.)

• Another idea: solving triangular systems is super fast! $$\left[\begin{array}{ccc} r_{11}&r_{12} & r_{13}\\0&r_{22}&r_{23}\\0&0&r_{33} \end{array} \right] \left[ \begin{array}{c}x_1\\x_2\\x_3\end{array} \right] = \left[ \begin{array}{c}b_1\\b_2\\b_3\end{array} \right]$$ Notice you immediately have $$x_{3} = b_3 / r_{33}$$. Then substitute $$x_3$$ and second row $$x_2r_{22}+x_3r_{23}=b_2$$ immediately gives you $$x_2$$. Computers can solve triangular systems at lightning speed applying this back substitutionbacksolve algorithm. (Don't reinvent the wheel, call Linear Algebra libraries to do this.) This is what Gordon is getting at when advocating taking a Cholesky decomposition on $$A$$. Once you have the Cholesky Decomposition $$A=RR'$$ where $$R$$ is lower triangular, you can solve $$A\mathbf{x}=\mathbf{b}$$ with two back substitutions, which is nearly free.

### Practical coding way forward

Since there's already an answer in R, I'll give some code in Matlab. In Matlab, there's the operator \ which basically examines your matrix and does the right thing to solve a linear system fast.

C_inv_X = C \ X;                     % Solves C^{-1}X
bhat = (X'*C_inv_X) \ (C_inv_X'*y);  % Final linear system solve


You could also:

R       = chol(C);  %Careful... does your library give you an UT or LT?
C_inv_X = R\(R'\X); %mldivide auto-detects triangular system
% in other environments, you may have to declare it to
% get full speed
bhat    = (X'*C_inv_X) \ (C_inv_X'*y);  % Final linear system solve


And then you have the triangular matrix $$R$$ around to use later.

### Big idea: solve linear systems rather than compute matrix inverse and multiply

Computing the inverse of a large matrix is slow and imprecise. If you have a huge matrix, don't do it. Instead:

• Observe $$A^{-1} \mathbf{b}$$ is the solution to $$A\mathbf{x} = \mathbf{b}$$. More generally, $$A^{-1}B$$ is the solution to $$AX=B$$. Anytime you see, $$A^{-1}B$$ you should think, "Hey, maybe I should solve linear systems rather than compute the matrix inverse?" (Note: for small problems, it probably doesn't matter.)

• Another idea: solving triangular systems is super fast! $$\left[\begin{array}{ccc} r_{11}&r_{12} & r_{13}\\0&r_{22}&r_{23}\\0&0&r_{33} \end{array} \right] \left[ \begin{array}{c}x_1\\x_2\\x_3\end{array} \right] = \left[ \begin{array}{c}b_1\\b_2\\b_3\end{array} \right]$$ Notice you immediately have $$x_{3} = b_3 / r_{33}$$. Then substitute $$x_3$$ and second row $$x_2r_{22}+x_3r_{23}=b_2$$ immediately gives you $$x_2$$. Computers can solve triangular systems at lightning speed applying this back substitution algorithm. (Don't reinvent the wheel, call Linear Algebra libraries to do this.) This is what Gordon is getting at when advocating taking a Cholesky decomposition on $$A$$. Once you have the Cholesky Decomposition $$A=RR'$$ where $$R$$ is lower triangular, you can solve $$A\mathbf{x}=\mathbf{b}$$ with two back substitutions, which is nearly free.

### Practical coding way forward

Since there's already an answer in R, I'll give some code in Matlab. In Matlab, there's the operator \ which basically examines your matrix and does the right thing to solve a linear system fast.

C_inv_X = C \ X;                     % Solves C^{-1}X
bhat = (X'*C_inv_X) \ (C_inv_X'*y);  % Final linear system solve


You could also:

R       = chol(C);  %Careful... does your library give you an UT or LT?
C_inv_X = R\(R'\X);
bhat    = (X'*C_inv_X) \ (C_inv_X'*y);  % Final linear system solve


And then you have the triangular matrix $$R$$ around to use later.

### Big idea: solve linear systems rather than compute matrix inverse and multiply

Computing the inverse of a large matrix is slow and imprecise. If you have a huge matrix, don't compute the inverse. Instead:

• Observe $$A^{-1} \mathbf{b}$$ is the solution to $$A\mathbf{x} = \mathbf{b}$$. More generally, $$A^{-1}B$$ is the solution to $$AX=B$$. Anytime you see, $$A^{-1}B$$ you should think, "Hey, maybe I should solve linear systems rather than compute the matrix inverse?" (Note: for small problems, it probably doesn't matter.)

• Another idea: solving triangular systems is super fast! $$\left[\begin{array}{ccc} r_{11}&r_{12} & r_{13}\\0&r_{22}&r_{23}\\0&0&r_{33} \end{array} \right] \left[ \begin{array}{c}x_1\\x_2\\x_3\end{array} \right] = \left[ \begin{array}{c}b_1\\b_2\\b_3\end{array} \right]$$ Notice you immediately have $$x_{3} = b_3 / r_{33}$$. Then substitute $$x_3$$ and second row $$x_2r_{22}+x_3r_{23}=b_2$$ immediately gives you $$x_2$$. Computers can solve triangular systems at lightning speed applying this backsolve algorithm. (Don't reinvent the wheel, call Linear Algebra libraries to do this.) This is what Gordon is getting at when advocating taking a Cholesky decomposition on $$A$$. Once you have the Cholesky Decomposition $$A=RR'$$ where $$R$$ is lower triangular, you can solve $$A\mathbf{x}=\mathbf{b}$$ with two back substitutions, which is nearly free.

### Practical coding way forward

Since there's already an answer in R, I'll give some code in Matlab. In Matlab, there's the operator \ which basically examines your matrix and does the right thing to solve a linear system fast.

C_inv_X = C \ X;                     % Solves C^{-1}X
bhat = (X'*C_inv_X) \ (C_inv_X'*y);  % Final linear system solve


You could also:

R       = chol(C);  %Careful... does your library give you an UT or LT?
C_inv_X = R\(R'\X); %mldivide auto-detects triangular system
% in other environments, you may have to declare it to
% get full speed
bhat    = (X'*C_inv_X) \ (C_inv_X'*y);  % Final linear system solve


And then you have the triangular matrix $$R$$ around to use later.

2 added 122 characters in body

### Big idea: solve linear systems rather than compute matrix inverse and multiply

Computing the inverse of a large matrix is slow and imprecise. If you have a huge matrix, don't do it. Instead:

• Observe $$A^{-1} \mathbf{b}$$ is the solution to $$A\mathbf{x} = \mathbf{b}$$. More generally, $$A^{-1}B$$ is the solution to $$AX=B$$. Anytime you see, $$A^{-1}B$$ you should think, "Hey, maybe I should solve a linear systemsystems rather than compute the matrix inverse?" (Note: for small problems, it probably doesn't matter.)

• Another idea: solving triangular systems is super fast! $$\left[\begin{array}{ccc} r_{11}&r_{12} & r_{13}\\0&r_{22}&r_{23}\\0&0&r_{33} \end{array} \right] \left[ \begin{array}{c}x_1\\x_2\\x_3\end{array} \right] = \left[ \begin{array}{c}b_1\\b_2\\b_3\end{array} \right]$$ Notice you immediately have $$x_{3} = b_3 / r_{33}$$. Then substitute $$x_3$$ and abovesecond row $$x_2r_{22}+x_3r_{23}=b_2$$ immediately gives you $$x_2$$. Computers can solve triangular systems at lightning speed applying this back substitution algorithm. (Don't reinvent the wheel, call Linear Algebra libraries to do this.) This is what Gordon is getting at when advocating taking a Cholesky decomposition on $$A$$. Once you have the Cholesky Decomposition $$A=RR'$$ where $$R$$ is lower triangular, you can solve $$A\mathbf{x}=\mathbf{b}$$ with two back substitutions, which is nearly free.

### Practical coding way forward

Since there's already an answer in R, I'll give some code in Matlab. In Matlab, there's the operator \ which basically examines your matrix and does the right thing to solve a linear system fast.

C_inv_X = C \ X;                     % Solves C^{-1}X
bhat = (X'*C_inv_X) \ (C_inv_X'*y);  % Final linear system solve


You could also:

R       = chol(C);  %Careful... does your library give you an UT or LT?
C_inv_X = R\(R'\X);
bhat    = (X'*C_inv_X) \ (C_inv_X'*y);  % Final linear system solve


And then you have the triangular matrix $$R$$ around to use later.

### Big idea: solve linear systems rather than compute matrix inverse and multiply

• $$A^{-1} \mathbf{b}$$ is the solution to $$A\mathbf{x} = \mathbf{b}$$. More generally, $$A^{-1}B$$ is the solution to $$AX=B$$. Anytime you see, $$A^{-1}B$$ you should think, "Hey, maybe I should solve a linear system rather than compute the matrix inverse?" (Note: for small problems, it probably doesn't matter.)

• Another idea: solving triangular systems is super fast! $$\left[\begin{array}{ccc} r_{11}&r_{12} & r_{13}\\0&r_{22}&r_{23}\\0&0&r_{33} \end{array} \right] \left[ \begin{array}{c}x_1\\x_2\\x_3\end{array} \right] = \left[ \begin{array}{c}b_1\\b_2\\b_3\end{array} \right]$$ Notice you immediately have $$x_{3} = b_3 / r_{33}$$. Then substitute $$x_3$$ and above row gives you $$x_2$$. Computers can solve triangular systems at lightning speed applying this back substitution algorithm. (Don't reinvent the wheel, call Linear Algebra libraries to do this.) This is what Gordon is getting at when advocating taking a Cholesky decomposition on $$A$$. Once you have the Cholesky Decomposition $$A=RR'$$ where $$R$$ is lower triangular, you can solve $$A\mathbf{x}=\mathbf{b}$$ with two back substitutions, which is nearly free.

### Practical coding way forward

In Matlab, there's the operator \ which basically examines your matrix and does the right thing to solve a linear system fast.

C_inv_X = C \ X;                     % Solves C^{-1}X
bhat = (X'*C_inv_X) \ (C_inv_X'*y);  % Final linear system solve


You could also:

R       = chol(C);  %Careful... does your library give you an UT or LT?
C_inv_X = R\(R'\X);
bhat    = (X'*C_inv_X) \ (C_inv_X'*y);  % Final linear system solve


And then you have the triangular matrix $$R$$ around to use later.

### Big idea: solve linear systems rather than compute matrix inverse and multiply

Computing the inverse of a large matrix is slow and imprecise. If you have a huge matrix, don't do it. Instead:

• Observe $$A^{-1} \mathbf{b}$$ is the solution to $$A\mathbf{x} = \mathbf{b}$$. More generally, $$A^{-1}B$$ is the solution to $$AX=B$$. Anytime you see, $$A^{-1}B$$ you should think, "Hey, maybe I should solve linear systems rather than compute the matrix inverse?" (Note: for small problems, it probably doesn't matter.)

• Another idea: solving triangular systems is super fast! $$\left[\begin{array}{ccc} r_{11}&r_{12} & r_{13}\\0&r_{22}&r_{23}\\0&0&r_{33} \end{array} \right] \left[ \begin{array}{c}x_1\\x_2\\x_3\end{array} \right] = \left[ \begin{array}{c}b_1\\b_2\\b_3\end{array} \right]$$ Notice you immediately have $$x_{3} = b_3 / r_{33}$$. Then substitute $$x_3$$ and second row $$x_2r_{22}+x_3r_{23}=b_2$$ immediately gives you $$x_2$$. Computers can solve triangular systems at lightning speed applying this back substitution algorithm. (Don't reinvent the wheel, call Linear Algebra libraries to do this.) This is what Gordon is getting at when advocating taking a Cholesky decomposition on $$A$$. Once you have the Cholesky Decomposition $$A=RR'$$ where $$R$$ is lower triangular, you can solve $$A\mathbf{x}=\mathbf{b}$$ with two back substitutions, which is nearly free.

### Practical coding way forward

Since there's already an answer in R, I'll give some code in Matlab. In Matlab, there's the operator \ which basically examines your matrix and does the right thing to solve a linear system fast.

C_inv_X = C \ X;                     % Solves C^{-1}X
bhat = (X'*C_inv_X) \ (C_inv_X'*y);  % Final linear system solve


You could also:

R       = chol(C);  %Careful... does your library give you an UT or LT?
C_inv_X = R\(R'\X);
bhat    = (X'*C_inv_X) \ (C_inv_X'*y);  % Final linear system solve


And then you have the triangular matrix $$R$$ around to use later.

1