Matrix Methods for

Computational Modeling
and Data Analytics

Virginia Tech · Spring 2019

Mark Embree
embree@[Link]

Chapter 2
Linear Systems from Resistor Networks
Ax=b
version of 21 January 2019

Circuits provide an elegant source of linear algebraic systems.

Here we shall only consider simple systems involving a battery and
some resistors, systems that sit at static equilibrium. As we proceed,
the goal is not so much a model of a given circuit, but a modeling
methodology that will apply to a broader set of problems we shall
study over the next few weeks.
These notes draw heavily, in spirit, details, and examples, from the
texts of Gilbert Strang1 and Steve Cox2 , and the lab experiments of 1
Gilbert Strang. Introduction to Applied
Cox et al.3 Mathematics. Wellesley-Cambridge
Press, Wellesley, MA, 1986
2
Steven J. Cox. Matrix Analysis in Situ.
Rice University, 2013
2.1 Resistor network modeling 3
Steven J. Cox, Mark Embree, and
Jeffrey M. Hokanson. Physical Laboratory
We begin with the example shown in Figure 2.1, consisting of six Manual for CAAM 335. Rice University,
2013
resistors and a constant voltage source. The goal is to determine what
the potential (voltage) is at three nodes in the network, x1 , x2 , x3 .

R1 x1 R3 x2 R5 x3 Figure 2.1: A circuit with an input

voltage and six resistors. We seek the
potential values at the nodes x1 , x2 , x3 .
Cox uses this circuit as a primitive
model for a neuron, where the horizon-
v0 R2 R4 R6 tal resistors model resistance caused
by the intercellular material, and the
vertical resistors model leakage through
the cell membrane.

While this network is simple, the methodology we shall derive ap-

plies to far more complicated circuits. By mastering this systematic
approach, you will develop skills of broad applicability – and the
linear algebra you need to solve such systems will prove even more
useful.
 11

With the constant voltage source v0 , this network sits at equi-

librium. The potential at x1 , x2 , and x3 will depend on v0 and the
relative strength of the resistors. We shall determine the decay of the
potential at points farther the voltage source in key modeling four
steps.

step 1 Compute voltage drops across resistors.

Across each of the six resistors, we compute the drop in voltage,
denoted e1 , . . . , e6 . As we consider the current flowing forth from
the voltage source v0 , we measure the voltage drop across R1 by
the potential before R1 minus the potential after, i.e.,

e1 = v 0 x1 .

We follow the same approach for the other five resistors. Since the
far side of R2 connects to ground,

e2 = x 1 0.

Similarly, the drops across R3 , . . . , R6 are

e3 = x 1 x2
e4 = x 2 0
e5 = x 2 x3
e6 = x 3 0.

Even for this simple network, these equations start to get tedious.
Just as for the population model in the last lecture, the organiza-
tion of individual equations into matrix-vector form will illumi-
nate. In this case we seek to relate six potential drops e1 , . . . , e6 to
three potential values x1 , x2 , and x3 , and we must handle the input
voltage as well. Collecting like terms in vectors, we have
2 3 2 3 2 32 3
x1
e1 v0 1 0 0
6 7 4 5
6 e2 7 6 0 7 6 1
6 7 6 7 6 0 0 7 x2
6e 7 6 0 7 6 1 7 x3
6 37 6 7 6 1 0 7
7
6 7=6 7 6 (2.1)
6 e4 7 6 0 7 6 0 1 0 7
7
6 7 6 7 6 7
4 e5 5 4 0 5 4 0 1 1 5
e6 0 0 0 1

which we denote as
e=v Ax. (2.2)

step 2 Apply Ohm’s Law.

We use Ohm’s Law to relate the voltage drop across each resistor
to current. You probably remember “V = IR” from physics class.
12

In this case we know “V” (the voltage drop) and “R” (the value of
the resistor), and seek “I,” so we use “I = V/R”. We shall write
the current at the six resistors as y1 , . . . , y6 . Then at each of the
resistors Ohm’s Law gives

y j = e j /R j , j = 1, . . . , 6.

As in Step 1, we want to write this in matrix–vector form, giving

2 3 2 32 3
y1 1/R1 0 0 0 0 0 e1
6 y2 7 6 0 1/R2 0 0 0 0 7 6 e2 7
6 7 6 76 7
6y 7 6 0 0 1/R3 0 0 0 7 6 7
6 37 6 7 6 e3 7
6 7=6 76 7,
6 y4 7 6 0 0 0 1/R4 0 0 7 6 e4 7
6 7 6 76 7
4 y5 5 4 0 0 0 0 1/R5 0 5 4 e5 5
y6 0 0 0 0 0 1/R6 e6

which we shall denote as

y = Ke. (2.3)

The direction of these currents is illustrated in Figure 2.2. The di-

agonal form of the matrix K corresponds to the fact that this step
of the modeling process does not encode any information about
the connectivity of the network: K just describes the material
properties of individual resistors.

R1 x1 R3 x2 R5 x3 Figure 2.2: The circuit from Figure 2.1,

with the current directions noted.
y1 y3 y5

v0 y2 R2 y4 R4 y6 R6

step 3 Apply Kirchhoff’s Current Law.

Having related the potential values to voltage drops, and voltage
drops to currents, we can now invoke the equilibrium condition
that will allow us to compute the unknown voltages at each node.
Kirchhoff’s Current Law says that the current entering each node
x1 , x2 , and x3 must sum to zero:

at x1 , y1 y2 y3 = 0;
at x2 , y3 y4 y5 = 0;
at x3 , y5 y6 = 0.
 13

Of course, we write this too as a matrix–vector product:

2 32 3 2 3
1 1 1 0 0 0 y1 0
6 7 6 y2 7 = 4 0 5.
4 0 0 1 1 1 0 56 7
6y 7 0
0 0 0 0 1 1 6 6 7
37
(2.4)
6 y4 7
6 7
4 y5 5
y6
Before giving this matrix a name, please pause to make this key
observation: both the voltage drop computation in Step 1 and
the current computation in Step 3 are determined by the wiring
paths in the circuit – they encode the structure of the circuit. The
first case maps the potential values x1 , x2 , and x3 (via Ohm’s law)
to the currents y1 , . . . , y6 . The second case does the reverse, in a
sense: it imposes a condition on the currents at each of the poten-
tials. Then it is no surprise then that the matrix in (2.4) is precisely
the transpose of the matrix A in (2.1). Thus we conserve notation by
writing (2.4) as
AT y = 0. (2.5)

step 4 Assembly.
Remember what we are after: given the voltage v0 (i.e., the vec-
tor v in (2.2), find the potentials x1 , x2 , and x3 . To obtain a clean
expression for these potential values, we need to assemble the
results of our first three steps.
Insert equation (2.3) for y into (2.5) to obtain

0 = AT y
= AT Ke.

Now insert equation (2.2) for e into this last result to obtain

0 = AT Ke
= AT K(v Ax).

Rewrite this equation, defining b := AT Kv, to get the the funda-

mental form:

AT KA x = b. (2.6)

Assuming we know values for the resistances R1 , . . . , R6 , we can

assemble the matrix AT KAx = b, and arrive at a simple linear
system of equations for the unknowns x1 , x2 , and x3 .
What size is the matrix AT KA? We see from the dimensions of the
ingredients ⇣ ⌘⇣ ⌘⇣ ⌘ ⇣ ⌘
3⇥6 6⇥6 6⇥3 = 3⇥3
14

that AT KA is a 3 ⇥ 3 matrix. Let us compute it, keeping symbolic

values for the resistances. First note that
2 32 3
1 1 1 0 0 0 1/R1 0 0 0 0 0
6 76 7
AT K = 4 0 0 1 1 1 0 56 0 1/R2 0 0 0 0 7
6 7
0 0 0 0 1 1 6 6 0 0 1/R3 0 0 0 7
7
6 0 0 0 1/R 0 0 7
6 4 7
6 7
4 0 0 0 0 1/R5 0 5
0 0 0 0 0 1/R6
2 3
1/R1 1/R2 1/R3 0 0 0
6 7
=4 0 0 1/R3 1/R4 1/R5 0 5.
0 0 0 0 1/R5 1/R6

Note: postmultiplying AT by the diagonal matrix K scaled the Yes, premultiplying by a diagonal matrix
scales the rows. Had we first computed
columns of AT by the diagonal entries. This is a general rule.
KA, we would have seen this.
Now compute (AT K)A:
2 3
1 0 0
2 36
6 1 0 0
7
7
1/R1 1/R2 1/R3 0 0 0 6 7
T 6 76 1 1 0 7
(A K)A = 4 0 0 1/R3 1/R4 1/R5 0 56
6
7
7
0 1 0
0 0 0 0 1/R5 1/R6 66
7
7
4 0 1 1 5
0 0 1
2 3
1/R1 + 1/R2 + 1/R3 1/R3 0
6 7
=6
4 1/R3 1/R3 + 1/R4 + 1/R5 1/R5 7.
5
0 1/R5 1/R5 + 1/R6

It remains to evaluate the right-hand side vector:

2 32 3 2 3
1/R1 1/R2 1/R3 0 0 0 v0 v0 /R1
6 76 7 6 7
b = (AT K)v = 4 0 0 1/R3 1/R4 1/R5 0 56 0 7 = 4 0 5.
6 7
0 0 0 0 1/R5 6
1/R6 6 0 7 7 0
6 0 7
6 7
6 7
4 0 5
0

We obtain the element-by-element version of the key equation (2.6):

2 32 3 2 3
1/R1 + 1/R2 + 1/R3 1/R3 0 x1 v0 /R1
6 76 7 6 7
6 1/R3 1/R3 + 1/R4 + 1/R5 1/R5 7 6 x2 7 = 6 0 7.
4 54 5 4 5
0 1/R5 1/R5 + 1/R6 x3 0

Notice that the matrix in this equation is symmetric. This must be the case, since

(AT KA) T = AT K T (AT ) T = AT KA.

(We used the fact that K = K T since K
is diagonal, and (AT ) T = A for all A.)
 15

What if all resistors have the same value, R Ohms? Then, clear-
ing 1/R terms, our equation becomes
2 32 3 2 3
3 1 0 x1 v0
4 1 3 1 5 4 x2 = 0 5 .
5 4 (2.7)
0 1 2 x3 0
Given a value for v0 , you can readily solve for x1 , x2 , and x3 .

student experiments

2.1. Consider the case (2.7) when all resistors have the same values,
R. Since this equation does not involve R, there is no way for the
resistors to influence the values of x1 , x2 , and x3 . Does this make
physical sense?

2.2. Implement the simple circuit in Figure 2.1 on a breadboard us-

ing a 9 V battery (or power supply, etc.) as your voltage source.
Use the same type of resistors for R1 , . . . , R6 , say R = 100 W. Mea-
sure the potentials x1 , x2 , and x3 at the nodes. Now solve (2.7) for
x1 , x2 , and x3 (in matlab using the \ command, or by hand using For example, if v0 = 9, use:
Gaussian elimination, as in the section below). How does your an- S = [3 -1 0; -1 3 -1 ; 0 -1 2];
b = [9; 0; 0];
swer compare to the prediction from our model? Now replace the x = S\b
vertical resistors with a different type of resistor, and repeat the
experiment.

2.3. Imagine the pattern established in the circuit of Figure 2.1 is

extended in a regular fashion to have 2N resistors (all of the same
value, R) with N 1 nodes x1 , . . . , x N . For example, N = 3 in
Figure 2.1, and Figure 2.3 shows the circuit for N = 5. How does
the N = 3 equation in (2.7) generalize to arbitrary larger values of
N? Can you deduce a general form for the matrix? If v0 and R are
held constant, what happens to the value of x N as N ! •? Does
this agree with your physical intuition?

R x1 R x2 R x3 R x4 R x5

v0 R R R R R

Figure 2.3: Extending the pattern of the

N = 3 circuit in Figure 2.1 to N = 5
2.4. Use the four steps developed above to build the linear system nodes and 2N resistors.

AT KAx = b for the branched circuit in Figure 2.4 (which models,

16

in a primitive but suggestive manner, a branched neuron.) Notice

how the branch affects the sparsity of the resulting matrix AT KA, The sparsity denotes the zero–nonzero
compared to the unbranched model for large N studied in Experi- structure of the matrix.

ment 2.3.

Figure 2.4: A circuit with an input

voltage and sixteen resistors in a
branched configuration. A neuron
would have many such branches
R6 R8 R10
modeling its dendrites.

R7 R9
R5 x3 x4 x5
R1 x1 R3 x2
R11

x6 R13 x7 R15 x8
v0 R2 R4

R12 R14 R16

2.2 Row reduction

Now it is time to solve our system for x1 , x2 , and x3 . For a concrete

example we shall address equation (2.7). The traditional way of solv-
ing such a system is Gaussian elimination. We presume you are al- “Gaussian” elimination is far more
ready familiar with this technique, but we will briefly recap it here. ancient than Gauss; an early user was
the 3rd century Chinese mathematician
Conventionally one writes the matrix and right-hand side together as Liu Hui. For more on this history, see:
the augmented matrix Joseph Grcar. How ordinary elimi-
nation became Gaussian elimination.
2 3
3 1 0 v0 Historia Math., 38:163–218, 2011
4 1 3 1 0 5
0 1 2 0

Standard Gaussian elimination proceeds by converting to zero the en-

tries below the diagonal of the matrix in the left of augmented form
through use of elementary row operations. These operations, which
are applied to entire rows of the augmented matrix, consist of three
techniques:

1. exchange two rows;

2. multiply a row by a nonzero scalar;

3. add one row to another row.

 17

These operations transform the augmented matrix while preserving Note that the apparently similar op-
the solution x1 , x2 , x3 (provided the operations are applied to both erations of exchanging, scaling, and
adding columns do not preserve x1 , x2 ,
the matrix and the right-hand side in the augmented form). We and x3 . Why?
demonstrate this technique on our 3 ⇥ 3 matrix, and, perhaps unlike
your past experience with Gaussian elimination, we shall keep the
variable term v0 in the right hand side of the equation.
Multiply row 2 by 3:
2 3
3 1 0 v0
4 3 9 3 0 5.
0 1 2 0

Add row 1 to row 2 to zero out the (2,1) entry:

2 3
3 1 0 v0
40 8 3 v0 5 .
0 1 2 0

Multiply row 3 by 8:
2 3
3 1 0 v0
40 8 3 v0 5 .
0 8 16 0

Add row 2 to row 3 to zero out the (3,2) entry:

2 3
3 1 0 v0
40 8 3 v0 5 . (2.8)
0 0 13 v0

Now the subdiagonal entries, i.e., those in the (2,1), (3,1), and (3,2)
positions, have been transformed to zero. The last augmented matrix
is equivalent to the linear system
2 32 3 2 3
3 1 0 x1 v0
40 8 3 5 4 x2 = v0 5 ,
5 4
0 0 13 x3 v0

and the upper triangular form of the matrix means that we can solve
the equations from the bottom–up, for the matrix–vector equation is
equivalent to the scalar equations

3x1 x2 + 0x3 = v0 (2.9)

8x2 3x3 = v0 (2.10)
13x3 = v0 (2.11)

First solve (2.11) for x3 :

v0
x3 = .
13
18

Substitute this formula for x3 into (2.10) and solve for x2 :

1⇣ 3⌘ 2v0
x2 = 1+ v0 = .
8 13 13

Finally, substitute the values of x1 and x2 into (2.9) and solve for x1 :

1⇣ 2⌘ 5v0
x1 = 1+ v0 = .
3 13 13

We summarize the solution:

2 3 2 3
x1 5
v 0 4 5
x = 4 x2 5 = 2 . (2.12)
13
x3 1

Now pause for one last essential step: Ask, Does the answer make The application gives this sanity check:
sense? First off, the units are correct: x1 , x2 , and x3 indeed have the made-up problems lack this safeguard.

same units as the voltage v0 . Moreover, with v0 > 0 all entries of x

are positive: we do not get negative voltages. And x1 > x2 > x3 : the
potentials decrease with distance from the voltage source. Increasing
v0 uniformly scales the potentials. All this agrees with our physical
intuition; the answer seems reasonable.

student experiments

2.5. Show how each of the three elementary row operations can be
encoded in the form of a matrix–matrix product.

(i) Design a matrix P j,k such that P j,k S swaps rows j and k of S.
For example, we want
2 3 2 3
a b c d e f
P1,2 4 d e f 5=4a b c 5.
g h i g h i

(ii) Design a matrix M j such that M j S multiplies row j of S by the

scalar g:
2 3 2 3
a b c a b c
M2 4 d e f 5 = 4 gd ge g f 5 .
g h i g h i

(iii) Design a matrix R j,k such that R j,k S replaces row j of S with
the sum of rows j and k:
2 3 2 3
a b c a b c
4
R3,1 d e 5
f = 4 d e f 5.
g h i a+g b+h c+i
 19

2.3 Gaussian elimination is LU factorization

Now is the time to develop a more mature understanding of Gaus-

sian elimination. This viewpoint was first articulated in the late
1940s, a time when early computer scientists were designing the first
codes to solve linear systems on computers.4 Focus on the matrix 4
G. W. Stewart. The decompositional
approach to matrix computation.
2 3 Computing Sci. Eng., 2:50–59, 2000
3 1 0
S= 4 1 3 15.
0 1 2

Since the (3,1) entry is already zero, our first task is to zero out the
(2,1) position. Varying slightly from the operations described above, Above we added the first row to three
we do so by adding 1/3 times the first row to the second row, consol- times the second row to preserve inte-
ger entries. When designing algorithms
idating two elementary row operations. Following on from Experi- for a computer, one prefers regularity to
ment 2.5, we perform this operation using a matrix-matrix product: cosmetic beauty. So, we always replace
a row by the sum of that row plus a
2 32 3 2 3 multiple of a different row.
1 0 0 3 1 0 3 1 0
4
L2,1 S = 1/3 1 0 5 4 1 3 5 4
1 = 0 8/3 15.
0 0 1 0 1 2 0 1 2

Next, manipulate L2,1 S to insert a zero in the (3,2) position, adding

3/8 of the new second row to the third row:
2 32 3 2 3
1 0 0 3 1 0 3 1 0
L3,2 (L2,1 S) = 4 0 1 0 5 4 0 8/3 1 5 = 4 0 8/3 1 5.
0 3/8 1 0 1 2 0 0 13/8

The resulting matrix is upper triangular (zero below the main diago-
nal), so we call it
2 3
3 1 0
U = 4 0 8/3 1 5.
0 0 13/8
Now since L3,2 L2,1 S = U, we can write

S = L2,11 L3,21 U. (2.13)

The inverses of the lower triangular matrices L3,2 and L2,1 are incredi-
bly simple:
2 3 2 3
1 0 0 1 0 0
L2,11 = 4 1/3 1 0 5 , L3,21 = 40 1 05. Check that the proposed inverses give
L2,1 L2,11 = I and L3,2 L3,21 = I.
0 0 1 0 3/8 1

Now compute the product L2,11 L3,21 in (2.13):

2 32 3 2 3
1 0 0 1 0 0 1 0 0
L2,11 L3,21 = 4 1/3 1 0 5 4 0 1 05 = 4 1/3 1 05.
0 0 1 0 3/8 1 0 3/8 1
20

This product of lower triangular matrices is also lower triangular, so

we call it 2 3
1 0 0
L := L2,11 L3,21 = 4 1/3 1 05. (2.14)
0 3/8 1

We arrive at S = L2,11 L3,21 U = LU, so

2 32 3 2 3
1 0 0 3 1 0 3 1 0
S = LU = 4 1/3 1 0 5 4 0 8/3 1 5 = 4 1 3 15.
0 3/8 1 0 0 13/8 0 1 2

This formula S = LU opens innumerable doors. For example, using

the properties of inverses described in Lecture 1,

1 1 1 1
S = (LU) =U L .

Maybe you are dubious that we have helped the situation by trans-
forming the problem of inverting one matrix S into the problem of
inverting two matrices U and L, then multiplying the results. But
since U and L are upper and lower triangular matrices, they are easy
to invert. For U we find
2 3
1/3 1/8 1/13
U 1=4 0 3/8 3/13 5
0 0 8/13

while for L, use equation (2.14) to see

2 3
1 0 0
1 Indeed, the inverses of upper/lower tri-
L = (L2,11 L3,21 ) 1
= L3,2 L2,1 = 4 1/3 1 05,
angular matrices are also upper/lower
1/8 3/8 1 triangular.

so we arrive at
2 32 3 2 3
1/3 1/8 1/13 1 0 0 5 2 1
1 4 5 4 5 1 4
S = 0 3/8 3/13 1/3 1 0 = 2 6 35.
13
0 0 8/13 1/8 3/8 1 1 3 8

More importantly, we can now use the LU factorization of S to solve

the equation Sx = b, i.e. (2.6): for one thing, we could write
2 32 3 2 3
5 2 1 v0 5
1 4 v0 4 5
x = S 1b = 2 6 354 0 5 = 2
13 13
1 3 8 0 1

in agreement with (2.12). However, we can arrive at x without explic-

itly computing S 1 . Use S = LU to arrive at

LUx = b,
 21

or

1
Ux = L b
= L3,2 L2,1 b.

Take a moment to savor this equation. Why? Because (aside from the
different way we chose to scale the rows when we added them), this
last equation is equivalent to the “augmented matrix” result in (2.8)!
When you applied elementary row operations to the augmented ma-
trix (starting with [S|b]), you reduced S to the upper triangular form
U, and applied those same reducing transformations to b, giving
L3,2 L2,1 b. When you then found x by solving equations (2.9)–(2.9),
you were inverting the triangular matrix U, i.e.,

1
x=U (L3,2 L2,1 )b
1 1
=U L b
1
=S b.

student experiments

2.6. We seek a general formula for the inverse of a matrix that has
simple structure like L2,1 and L3,2 , i.e., an identity matrix with a
single off-diagonal entry set to a (the others being zero). We can
write such a matrix as
I + ae j ekT ,

where a is the entry that goes in the ( j, k ) position. (Recall from

Lecture 1 that e` denotes the `th column of the identity matrix, so
e j ekT is the matrix that is zero in all entries, save for a 1 in the (j,k)
position.) Inspired by the form of L2,11 and L3,21 , guess a formula
for (I + ae j ek ) 1 (for an arbitrary dimension) and show that it
works.

2.4 Epilogue: pivoted LU factorizations

Thinking of Gaussian elimination as the matrix factorization S = LU

is a higher form of thinking that has broad consequences for both
theory and numerical algorithms. However, not every matrix has a
decomposition of this form. For example, if

0 1
S= ,
1 1

then S is an invertible matrix, but there is no way to write S = LU

with L lower triangular and U upper triangular. If you tried to apply
22

the style of Gaussian elimination you learned in high school to this

equation, you would start by swapping the rows. That is the key
problem: row swaps destroy the lower-triangular structure in L. In a
course in numerical analysis, you will learn that you can encode the
row swaps by premultiplying S by a matrix P whose columns are the
same as the identity matrix, but arranged in a different order. For any
invertible matrix S we can always factor

PS = LU.
For more, see, e.g., the textbook:
There is so much more to say about this factorization, but in this Lloyd N. Trefethen and David Bau,
course we must move on now to other topics. . . . III. Numerical Linear Algebra. SIAM,
Philadelphia, 1997