Chapter 6 Experiment 6: Step-by-Step Elimination in Gaussian Elimination/PLU Decomposition and the Rank-One View
Goal of the experiment
Use an interactive tool to watch, step by step, how takes shape, and understand, from the point of view of the rank-one decomposition, the cumulative meaning of .
Part 1: Interactive Exploration¶
The Six Panels¶
| Position | Panel | Description |
|---|---|---|
| Top left | (original) | Unchanged; the reference throughout |
| Top middle | (permutation) | Records every row exchange; orange box = current pivot position |
| Top right | rank-one accumulation | Sum of the completed rank-one layers, gradually approaching |
| Bottom left | (lower triangular) | The elimination multipliers are filled in entry by entry |
| Bottom middle | (row echelon) | Reduced step by step to REF; teal dashed box = the row being eliminated |
| Bottom right | Current rank-one layer | The contributed in this round |
Color meaning: light blue = negative, white = zero, the panel’s own light warm color = positive.
How to Use the Tool¶
Click “▶ Decompose” to load the default matrix
Step forward with the ◀ / ▶ buttons or by dragging the slider
You can edit the values in the entry boxes, or click “Random” to explore different cases
“Resize” switches between and (non-square matrices included)
Source
# ══ Interactive PLU decomposition teaching widget (formerly plu_widget.py, now merged into this notebook) ══
"""
Interactive PLU decomposition teaching widget
Supports m×n matrices (3≤m≤5, 3≤n≤5; non-square matrices allowed)
Textbook: Linear Algebra with Python, §6.3 Gaussian elimination and the PLU decomposition
"""
import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots
import ipywidgets as widgets
from IPython.display import display, clear_output
import warnings
warnings.filterwarnings("ignore")
np.set_printoptions(precision=4, suppress=True, linewidth=100)
# ─── Color constants: book palette ───────────────────────────────────────────────
C_BG = "#F8F8F8"
C_GRID = "#D6D6D6"
C_AXIS = "#000000"
C_V1 = "#57068C"
C_V2 = "#006385"
C_T1 = "#2AD2C9"
C_T2 = "#8900E1"
C_WARN = "#FF5D47"
C_AUX = "#AB82C5"
RANK1_COLORS = [C_V1, C_WARN, C_V2, C_T1, C_T2, C_AUX]
# ─── PLU decomposition (with step-by-step snapshots) ───────────────────────────────────────────────
def plu_with_steps(A, tol=1e-12, skip_pivot=False):
# skip_pivot=True: no row exchanges (no partial pivoting); the current diagonal entry is the pivot.
# On a zero pivot with nonzero entries below it, record an explanatory step and stop.
m, n = A.shape
P = np.eye(m)
L = np.eye(m)
U = A.copy().astype(float)
steps = []
def snap(label, pivot=None, elim_row=None, rank1_idx=0):
steps.append(dict(label=label, P=P.copy(), L=L.copy(), U=U.copy(),
pivot=pivot, elim_row=elim_row, rank1_idx=rank1_idx))
snap("Initial state")
cur_row = cur_col = pivot_count = 0
while cur_row < m and cur_col < n:
col_seg = U[cur_row:, cur_col]
rel_best = 0 if skip_pivot else np.argmax(np.abs(col_seg))
best_row = cur_row + rel_best
if abs(U[best_row, cur_col]) < tol:
if skip_pivot and np.max(np.abs(col_seg)) >= tol:
snap(f"Step {pivot_count+1}b: pivot U[{cur_row},{cur_col}]=0; without pivoting the multiplier is ∞ (division by zero), so elimination cannot proceed and the decomposition stops",
pivot=(cur_row, cur_col), rank1_idx=pivot_count)
return steps, P, L, U
cur_col += 1
continue
if best_row != cur_row:
U[[cur_row, best_row]] = U[[best_row, cur_row]]
P[[cur_row, best_row]] = P[[best_row, cur_row]]
if cur_row > 0:
L[[cur_row, best_row], :cur_row] = L[[best_row, cur_row], :cur_row]
snap(f"Step {pivot_count+1}a: row exchange, row {cur_row} ↔ row {best_row}",
pivot=(cur_row, cur_col), rank1_idx=pivot_count)
pivot = U[cur_row, cur_col]
for i in range(cur_row + 1, m):
factor = U[i, cur_col] / pivot
if abs(factor) > tol:
L[i, cur_row] = factor
U[i, :] -= factor * U[cur_row, :]
snap(f"Step {pivot_count+1}b: eliminate row {i}, multiplier={factor:.4f}",
pivot=(cur_row, cur_col), elim_row=i, rank1_idx=pivot_count)
pivot_count += 1
snap(f"Round {pivot_count} of elimination done (pivot=({cur_row},{cur_col}))",
pivot=(cur_row, cur_col), rank1_idx=pivot_count)
cur_row += 1
cur_col += 1
snap("PLU decomposition complete", rank1_idx=pivot_count)
return steps, P, L, U
def compute_rank1_layers(L, U, P):
m, n = U.shape
rank = sum(1 for i in range(min(m, n)) if abs(U[i, i]) > 1e-12)
layers = []
cumsum = np.zeros((m, n))
for k in range(rank):
layer = L[:, k:k+1] @ U[k:k+1, :]
cumsum = cumsum + layer
layers.append(dict(rank1=layer.copy(), cumsum=(P.T @ cumsum).copy()))
return layers
# ─── Plotting ─────────────────────────────────────────────────────────────────
def build_figure(step_data, rank1_layers, A_orig, step_idx, n_steps):
m, n = A_orig.shape
P_snap = step_data["P"]
L_snap = step_data["L"]
U_snap = step_data["U"]
pivot = step_data["pivot"]
elim_row = step_data["elim_row"]
ri = step_data["rank1_idx"]
if ri > 0 and ri <= len(rank1_layers):
cum_mat = rank1_layers[ri-1]["cumsum"]
r1_mat = rank1_layers[ri-1]["rank1"]
else:
cum_mat = np.zeros((m, n))
r1_mat = np.zeros((m, n))
# Coloring logic: negative → light blue (shared), zero → white, positive → the panel's own light warm color
# Three-stop symmetric colorscale: stop 0 = negative color, stop 0.5 = white (zero), stop 1 = positive color
NEG = "#AECFE0" # shared color for negatives: light blue (meaning: negative)
WH = "#FFFFFF" # zero: white
CS_neutral = [[0, NEG], [0.5, WH], [1, "#D0B8E8"]] # A/P: light violet for positives
CS_L = [[0, NEG], [0.5, WH], [1, "#C8B0E8"]] # L: light violet for positives
CS_U = [[0, NEG], [0.5, WH], [1, "#F4C8A0"]] # U: light orange for positives (clearly distinct from the light blue of negatives)
_r1pos = ["#C8B0E8","#FFCCC8","#A8CDE0","#A8EAE8","#D0A8E8","#DCCCE8"]
CS_rank1 = [[0, NEG], [0.5, WH], [1, _r1pos[min(ri, 5)]]]
r1_title = (f"Rank-one layer {ri}: c<sub>{ri-1}</sub>·r<sub>{ri-1}</sub><sup>T</sup>"
if ri > 0 else "Rank-one layer (not started yet)")
fig = make_subplots(
rows=2, cols=3,
subplot_titles=[
"<b>A</b> (original)", "<b>P</b> (permutation)",
f"<b>P<sup>T</sup>LU</b> rank-one accumulation (first {ri} layers)",
"<b>L</b> (lower triangular)", "<b>U</b> (row echelon)", r1_title,
],
vertical_spacing=0.14, horizontal_spacing=0.06,
)
# Panel definitions: (matrix, colorscale, pivot position, row being eliminated, row, col)
panels = [
(A_orig, CS_neutral, None, None, 1, 1),
(P_snap, CS_neutral, pivot, None, 1, 2),
(cum_mat, CS_rank1, None, None, 1, 3),
(L_snap, CS_L, pivot, elim_row,2, 1),
(U_snap, CS_U, pivot, elim_row,2, 2),
(r1_mat, CS_rank1, None, None, 2, 3),
]
# Plotly subplot axis names: the first uses "x"/"y", the later ones "x2"/"y2"...
xref_map = {(1,1):"x", (1,2):"x2",(1,3):"x3",
(2,1):"x4",(2,2):"x5", (2,3):"x6"}
yref_map = {(1,1):"y", (1,2):"y2",(1,3):"y3",
(2,1):"y4",(2,2):"y5", (2,3):"y6"}
all_shapes = []
# P is an m×m square matrix and so is L; only A/U/cum/r1 are m×n
panel_shapes = {
(1,1): (m, n),
(1,2): (m, m),
(1,3): (m, n),
(2,1): (m, m),
(2,2): (m, n),
(2,3): (m, n),
}
for mat, cs, hpv, hr, row, col in panels:
pr, pc = panel_shapes[(row, col)]
zmax = max(np.max(np.abs(mat)), 1e-9)
fig.add_trace(go.Heatmap(
z=mat.tolist(), zmin=-zmax, zmax=zmax, zmid=0,
colorscale=cs, showscale=False,
hovertemplate="row %{y}, column %{x}<br>value=%{z:.4f}<extra></extra>",
), row=row, col=col)
xref = xref_map[(row, col)]
yref = yref_map[(row, col)]
# Value labels (appended with add_annotation, so the subplot_titles are not overwritten)
for i in range(pr):
for j in range(pc):
v = mat[i, j]
txt = f"{v:.2f}" if abs(v) >= 0.005 else "0"
fig.add_annotation(
x=j, y=i, text=txt, showarrow=False,
font=dict(color="#111111", size=10, family="monospace"),
xref=xref, yref=yref,
)
# Highlight box for the pivot
if hpv is not None:
pr2, pc2 = hpv
all_shapes.append(dict(
type="rect", x0=pc2-.5, x1=pc2+.5, y0=pr2-.5, y1=pr2+.5,
line=dict(color=C_WARN, width=3),
fillcolor="rgba(0,0,0,0)", xref=xref, yref=yref,
))
# Highlight box for the row being eliminated
if hr is not None:
all_shapes.append(dict(
type="rect", x0=-.5, x1=pc-.5,
y0=hr-.5, y1=hr+.5,
line=dict(color=C_T1, width=2, dash="dot"),
fillcolor="rgba(42,210,201,0.08)", xref=xref, yref=yref,
))
# Axis ticks
ax_kw = dict(showgrid=False, zeroline=False, tickfont=dict(size=8))
ax_updates = {}
for (row, col), xref in xref_map.items():
xax = "xaxis" if xref == "x" else f"xaxis{xref[1:]}"
yax = "yaxis" if yref_map[(row,col)] == "y" else f"yaxis{yref_map[(row,col)][1:]}"
pr, pc = panel_shapes[(row, col)]
ax_updates[xax] = dict(**ax_kw,
tickvals=list(range(pc)),
ticktext=[f"col {j}" for j in range(pc)])
ax_updates[yax] = dict(**ax_kw,
tickvals=list(range(pr)),
ticktext=[f"row {i}" for i in range(pr)],
autorange="reversed")
fig.update_layout(
height=max(420, 120 * m),
paper_bgcolor=C_BG, plot_bgcolor=C_BG,
margin=dict(l=10, r=10, t=80, b=10),
font=dict(family="Microsoft JhengHei,SimHei,sans-serif",
color=C_AXIS, size=11),
shapes=all_shapes,
showlegend=False,
title=dict(
text=(f"<b>PLU decomposition</b> "
f"<span style='color:{C_WARN}'>{step_data['label']}</span> "
f"<span style='color:{C_AUX};font-size:11px'>"
f"step {step_idx}/{n_steps-1}</span> "
f"<span style='color:{C_V2};font-size:11px'>"
f"matrix {m}×{n}</span>"),
x=0.5, xanchor="center",
font=dict(size=13, color=C_V1),
),
**ax_updates,
)
return fig
# ─── Widget ───────────────────────────────────────────────────────────────
class PLUWidget:
DEFAULT_A = np.array([
[2., 1., 1.],
[4., -6., 0.],
[-2., 7., 2.],
])
def __init__(self, A=None):
self.A = (A if A is not None else self.DEFAULT_A).copy().astype(float)
self._validate_shape()
self._build_steps()
self._build_ui()
def _validate_shape(self):
m, n = self.A.shape
m = max(3, min(5, m))
n = max(3, min(5, n))
if self.A.shape != (m, n):
new = np.zeros((m, n))
rm, rn = min(m, self.A.shape[0]), min(n, self.A.shape[1])
new[:rm, :rn] = self.A[:rm, :rn]
self.A = new
def _build_steps(self):
self.steps, self.P_f, self.L_f, self.U_f = plu_with_steps(self.A)
self.rank1_layers = compute_rank1_layers(self.L_f, self.U_f, self.P_f)
# ── Rebuild the entry boxes ───────────────────────────────────────────────────────
def _rebuild_grid(self, A):
"""Rebuild entry_widgets from the matrix A and update grid_box.children."""
m, n = A.shape
self.entry_widgets = []
rows = []
cell_w = "58px" if n >= 5 else "68px"
for i in range(m):
row_ws = []
for j in range(n):
w = widgets.FloatText(
value=A[i, j],
layout=widgets.Layout(width=cell_w),
style={"description_width": "0px"},
)
row_ws.append(w)
self.entry_widgets.append(w)
rows.append(widgets.HBox(row_ws))
self.grid_box.children = rows
def _size_label_html(self, m, n):
return (f"<span style='color:{C_V2};font-size:12px'>"
f"A: {m} rows × {n} columns (P, L: {m}×{m}; U: {m}×{n})</span>")
# ── Building the UI ──────────────────────────────────────────────────────────
def _build_ui(self):
m, n = self.A.shape
# Size selectors
self.sel_m = widgets.Dropdown(
options=[(f"{v} rows", v) for v in [3,4,5]], value=m,
layout=widgets.Layout(width="100px"),
style={"description_width": "0px"},
)
self.sel_n = widgets.Dropdown(
options=[(f"{v} cols", v) for v in [3,4,5]], value=n,
layout=widgets.Layout(width="100px"),
style={"description_width": "0px"},
)
btn_resize = widgets.Button(
description="Resize",
layout=widgets.Layout(width="90px"),
button_style="warning",
)
btn_resize.on_click(self._on_resize)
size_row = widgets.HBox(
[widgets.HTML("<span style='font-size:12px;line-height:28px'>Rows:</span>"),
self.sel_m,
widgets.HTML("<span style='font-size:12px;line-height:28px;margin-left:6px'>Columns:</span>"),
self.sel_n,
widgets.HTML("<span style='margin-left:6px'></span>"),
btn_resize],
layout=widgets.Layout(align_items="center", flex_wrap="wrap"),
)
self.size_label = widgets.HTML(value=self._size_label_html(m, n))
self.grid_box = widgets.VBox([])
self._rebuild_grid(self.A)
self.btn_default = widgets.Button(description="Reset",
layout=widgets.Layout(width="90px"))
self.btn_random = widgets.Button(description="Random",
layout=widgets.Layout(width="90px"))
self.btn_apply = widgets.Button(description="▶ Decompose",
button_style="success",
layout=widgets.Layout(width="110px"))
self.btn_default.on_click(self._on_reset)
self.btn_random.on_click(self._on_random)
self.btn_apply.on_click(self._on_apply)
input_panel = widgets.VBox([
widgets.HTML(f"<b style='color:{C_V1};font-size:13px'>Matrix size and entries</b>"),
size_row,
self.size_label,
self.grid_box,
widgets.HBox([self.btn_default, self.btn_random, self.btn_apply],
layout=widgets.Layout(margin="6px 0 0 0")),
], layout=widgets.Layout(
border=f"1.5px solid {C_AUX}", padding="10px",
border_radius="6px", width="330px",
))
# Slider
self.slider = widgets.IntSlider(
min=0, max=len(self.steps)-1, value=0,
description="Step",
continuous_update=False,
layout=widgets.Layout(width="90%"),
style={"description_width": "40px"},
)
self.btn_prev = widgets.Button(description="◀ Back",
layout=widgets.Layout(width="100px"))
self.btn_next = widgets.Button(description="Next ▶",
layout=widgets.Layout(width="100px"))
self.btn_prev.on_click(self._on_prev)
self.btn_next.on_click(self._on_next)
controls = widgets.HBox(
[self.btn_prev, self.slider, self.btn_next],
layout=widgets.Layout(align_items="center", width="100%"),
)
self.plot_out = widgets.Output()
self.rank1_html = widgets.HTML(value=self._rank1_html(0))
self.final_html = widgets.HTML(value=self._final_html())
self._redraw(0)
self.slider.observe(self._on_slider, names="value")
right_panel = widgets.VBox(
[controls, self.plot_out],
layout=widgets.Layout(flex="1"),
)
self.layout = widgets.VBox([
widgets.HTML(
f"<h3 style='color:{C_V1};margin:4px 0'>Interactive PLU Decomposition Tutorial</h3>"
f"<p style='color:#666;font-size:12px;margin:0'>"
f"Supports matrices with 3–5 rows × 3–5 columns (non-square included)</p>"
),
widgets.HBox(
[input_panel, right_panel],
layout=widgets.Layout(align_items="flex-start", gap="10px"),
),
self.rank1_html,
self.final_html,
])
# ── Events ─────────────────────────────────────────────────────────────
def _on_resize(self, _):
new_m = self.sel_m.value
new_n = self.sel_n.value
old_m, old_n = self.A.shape
# keep whatever values can be kept and fill the rest with 0
new_A = np.zeros((new_m, new_n))
rm = min(new_m, old_m)
rn = min(new_n, old_n)
new_A[:rm, :rn] = self.A[:rm, :rn]
self.A = new_A
self._rebuild_grid(self.A)
self.size_label.value = self._size_label_html(new_m, new_n)
def _on_prev(self, _):
if self.slider.value > 0:
self.slider.value -= 1
def _on_next(self, _):
if self.slider.value < self.slider.max:
self.slider.value += 1
def _on_slider(self, change):
idx = change["new"]
self._redraw(idx)
self.rank1_html.value = self._rank1_html(idx)
def _on_apply(self, _):
m, n = self.A.shape
new_A = np.array([
[self.entry_widgets[i*n + j].value for j in range(n)]
for i in range(m)
], dtype=float)
self.A = new_A
self._build_steps()
self.slider.max = len(self.steps) - 1
self.slider.value = 0
self.final_html.value = self._final_html()
self._redraw(0)
self.rank1_html.value = self._rank1_html(0)
def _on_reset(self, _):
m, n = self.A.shape
new_A = np.zeros((m, n))
dm, dn = self.DEFAULT_A.shape
rm, rn = min(m, dm), min(n, dn)
new_A[:rm, :rn] = self.DEFAULT_A[:rm, :rn]
self.A = new_A
self._rebuild_grid(self.A)
self._on_apply(None)
def _on_random(self, _):
m, n = self.A.shape
self.A = np.random.randint(-5, 6, size=(m, n)).astype(float)
self._rebuild_grid(self.A)
self._on_apply(None)
# ── Redrawing ─────────────────────────────────────────────────────────────
def _redraw(self, idx):
fig = build_figure(
self.steps[idx], self.rank1_layers,
self.A, idx, len(self.steps),
)
with self.plot_out:
clear_output(wait=True)
fig.show()
# ── Rank-one HTML ────────────────────────────────────────────────────────
def _rank1_html(self, idx):
ri = self.steps[idx]["rank1_idx"]
total = len(self.rank1_layers)
if total == 0:
return ""
pct = int(ri / total * 100) if total else 0
bar = (f"<div style='background:{C_GRID};border-radius:4px;"
f"height:8px;width:100%;margin:4px 0'>"
f"<div style='background:{C_V1};width:{pct}%;height:100%;"
f"border-radius:4px'></div></div>")
items = " ".join(
f"<span style='color:{RANK1_COLORS[k%len(RANK1_COLORS)]};"
f"opacity:{'1.0' if k<ri else '0.35'};margin-right:8px'>"
f"{'✓' if k<ri else ('▶' if k==ri-1 else '○')} "
f"layer {k}: c<sub>{k}</sub>r<sub>{k}</sub><sup>T</sup></span>"
for k in range(total)
)
parts = " + ".join(
f"<span style='color:{RANK1_COLORS[k%len(RANK1_COLORS)]}'>"
f"c<sub>{k}</sub>r<sub>{k}</sub><sup>T</sup></span>"
for k in range(total)
)
return (
f"<div style='border:1.5px solid {C_AUX};border-radius:6px;"
f"padding:10px 14px;background:{C_BG};font-size:13px'>"
f"<b style='color:{C_V1}'>Rank-one decomposition: </b>"
f"A = P<sup>T</sup>({parts}) "
f"completed {ri}/{total} layers"
f"{bar}<div>{items}</div></div>"
)
# ── Final-result HTML ────────────────────────────────────────────────────
def _final_html(self):
def mat2html(M, color):
rows = "".join(
"<tr>" + "".join(
f"<td style='padding:2px 6px;text-align:right;"
f"font-family:monospace;font-size:11px'>{v:.3f}</td>"
for v in row
) + "</tr>"
for row in M
)
return (f"<table style='border:1.5px solid {color};"
f"border-collapse:collapse'>{rows}</table>")
err = np.max(np.abs(self.P_f @ self.A - self.L_f @ self.U_f))
check = ("✓ PA=LU verified" if err < 1e-9
else f"✗ error={err:.2e}")
color = "green" if err < 1e-9 else "red"
m, n = self.A.shape
return (
f"<div style='border:1.5px solid {C_V1};border-radius:6px;"
f"padding:10px 14px;background:{C_BG};margin-top:6px;font-size:13px'>"
f"<b style='color:{C_V1}'>Final result PA = LU</b> "
f"<span style='color:{C_AUX};font-size:11px'>"
f"P:{m}×{m}, L:{m}×{m}, U:{m}×{n}</span><br>"
f"<div style='display:flex;align-items:center;gap:8px;"
f"margin-top:8px;flex-wrap:wrap'>"
f"<div><b>P</b><br>{mat2html(self.P_f, C_V2)}</div>"
f"<span style='font-size:18px;color:{C_AUX}'>×</span>"
f"<div><b>A</b><br>{mat2html(self.A, C_AXIS)}</div>"
f"<span style='font-size:18px;color:{C_AUX}'>=</span>"
f"<div><b>L</b><br>{mat2html(self.L_f, C_T2)}</div>"
f"<span style='font-size:18px;color:{C_AUX}'>×</span>"
f"<div><b>U</b><br>{mat2html(self.U_f, C_V2)}</div>"
f"<div style='margin-left:10px'>"
f"<span style='color:{color}'>{check}</span><br>"
f"<span style='color:#777;font-size:10px'>‖PA−LU‖∞={err:.2e}</span>"
f"</div></div></div>"
)
def show(self):
display(self.layout)
def plu_explorer(A=None):
"""
Launch the interactive PLU decomposition explorer.
Parameters
----------
A : ndarray, optional
Initial matrix (3≤m≤5, 3≤n≤5). Defaults to the 3×3 example matrix.
Examples
--------
>>> plu_explorer()
>>> plu_explorer(np.array([[1,2,3,4],[0,4,5,6],[1,8,0,2],[3,1,2,0]], dtype=float))
"""
w = PLUWidget(A)
w.show()
return w
# ══ Display settings for this experiment ══
np.set_printoptions(precision=4, suppress=True, linewidth=80)
Experiment A: The Standard Case (Row Exchanges Needed)¶
Press “Next” step by step and observe how row exchanges change and how the elimination multipliers are filled into .
A_demo = np.array([
[ 2., 1., 1.],
[ 4., -6., 0.],
[-2., 7., 2.],
])
w_A = plu_explorer(A_demo)
Loading...
Experiment B: A Rank-Deficient Matrix ()¶
Note that row 0 and row 1 are proportional. Observe the zero row of after elimination, and how the rank-one accumulation stops growing after a certain step.
B_demo = np.array([
[1., 2., 3.],
[2., 4., 6.],
[0., 1., 2.],
])
w_B = plu_explorer(B_demo)
Loading...
Experiment C: A Non-Square Matrix ()¶
Note that and are still ; only is . Observe how the number of rank-one layers is related to the number of pivots of .
C_demo = np.array([
[1., 2., 0., 3.],
[2., 4., 1., 7.],
[1., 1., 2., 4.],
])
w_C = plu_explorer(C_demo)
Loading...
Part 2: Observation and Reflection¶
After completing the interactive exploration of Part 1, answer the following questions based on your observations. We suggest that you first try to answer in your own words and then run the verification code to confirm.
Question 1: The Structure of the Permutation Matrix ¶
In Experiment A, what kind of matrix is the final ?
(a) Run the code below and observe how many nonzero entries each row and each column of has.
(b) Does satisfy ? Verify it with the code, and explain the geometric meaning of this property.
(c) Is the result of the PLU decomposition unique? If all the leading principal minors of the matrix are nonzero (this guarantees an LU decomposition without row exchanges; but this tool uses partial pivoting and may still exchange rows, as in Experiment A), and no row exchanges are made during the decomposition, what is ? To which decomposition does PLU then reduce?
# Question 1: verification
steps_A, P, L, U = plu_with_steps(A_demo)
print("P =")
print(P)
print()
# (b) Verify P^T P = I
PtP = P.T @ P
print("P^T P =")
print(PtP)
print(f"P^T P is the identity matrix: {np.allclose(PtP, np.eye(len(P)))}")
print()
# Number of nonzero entries in each row/column
print("Nonzero entries in each column of P:", [int(np.sum(P[:, j] != 0)) for j in range(P.shape[1])])
print("Nonzero entries in each row of P:", [int(np.sum(P[i, :] != 0)) for i in range(P.shape[0])])
P =
[[0. 1. 0.]
[1. 0. 0.]
[0. 0. 1.]]
P^T P =
[[1. 0. 0.]
[0. 1. 0.]
[0. 0. 1.]]
P^T P is the identity matrix: True
Nonzero entries in each column of P: [1, 1, 1]
Nonzero entries in each row of P: [1, 1, 1]
Your observations (fill in):
(a)
(b)
(c)
Question 2: The Diagonal Entries of ¶
(a) Run the code below and read off the diagonal entries of the final in Experiment A. What pattern do their values follow?
(b) The elimination multiplier is (the values before elimination step ). If the pivot is very small (close to zero), what goes wrong with ? This is exactly why “partial pivoting” chooses the entry of largest absolute value as the pivot; explain your understanding.
(c) For the rank-deficient matrix of Experiment B, are the diagonal entries of all equal to 1?
# Question 2: verification
print("L of Experiment A:")
print(L)
print(f"Diagonal entries of L: {np.diag(L)}")
print()
# Experiment B
steps_B, P_B, L_B, U_B = plu_with_steps(B_demo)
print("L of Experiment B (rank-deficient):")
print(L_B)
print(f"Diagonal entries of L_B: {np.diag(L_B)}")
print()
print("U of Experiment B:")
print(U_B)
print(f"Pivots of U (diagonal entries): {np.diag(U_B)}")
L of Experiment A:
[[ 1. 0. 0. ]
[ 0.5 1. 0. ]
[-0.5 1. 1. ]]
Diagonal entries of L: [1. 1. 1.]
L of Experiment B (rank-deficient):
[[1. 0. 0. ]
[0. 1. 0. ]
[0.5 0. 1. ]]
Diagonal entries of L_B: [1. 1. 1.]
U of Experiment B:
[[2. 4. 6.]
[0. 1. 2.]
[0. 0. 0.]]
Pivots of U (diagonal entries): [2. 1. 0.]
Your observations (fill in):
(a)
(b)
(c)
Question 3: The Rank-One View—How Much “Information” the Decomposition Carries¶
Each rank-one layer is a matrix of rank 1.
(a) In Experiment A, after how many layers does the “rank-one accumulation” panel at the top right become identical to the original matrix ? Run the code below to compare the error layer by layer.
(b) For Experiment B (), what does rank-one layer 2 (, if it exists) contribute? That is, how does the accumulated matrix change after layer 2 is added?
(c) Think about it: how is the number of rank-one layers related to the rank of the matrix? State your conjecture in terms of .
# Question 3: verification
layers_A = compute_rank1_layers(L, U, P)
print("Experiment A: error of the accumulation, layer by layer, ‖accumulation - A‖∞")
print("-" * 40)
for k, layer in enumerate(layers_A):
err = np.max(np.abs(layer["cumsum"] - A_demo))
print(f" error after accumulating the first {k+1} layers = {err:.6f}")
print()
# Experiment B
layers_B = compute_rank1_layers(L_B, U_B, P_B)
print(f"Experiment B (rank-deficient): number of rank-one layers = {len(layers_B)}")
print("Error of the accumulation, layer by layer, ‖accumulation - B‖∞")
print("-" * 40)
for k, layer in enumerate(layers_B):
err = np.max(np.abs(layer["cumsum"] - B_demo))
print(f" error after accumulating the first {k+1} layers = {err:.6f}")
# Non-square matrix
layers_C = compute_rank1_layers(*plu_with_steps(C_demo)[1:])
print(f"\nExperiment C (3×4): number of rank-one layers = {len(layers_C)}")
Experiment A: error of the accumulation, layer by layer, ‖accumulation - A‖∞
----------------------------------------
error after accumulating the first 1 layers = 4.000000
error after accumulating the first 2 layers = 1.000000
error after accumulating the first 3 layers = 0.000000
Experiment B (rank-deficient): number of rank-one layers = 2
Error of the accumulation, layer by layer, ‖accumulation - B‖∞
----------------------------------------
error after accumulating the first 1 layers = 2.000000
error after accumulating the first 2 layers = 0.000000
Experiment C (3×4): number of rank-one layers = 3
Your observations (fill in):
(a)
(b)
(c) My conjecture: the number of rank-one layers = ______
Question 4: The Effect of Row Exchanges on the Decomposition¶
Construct a matrix with a zero diagonal entry (but invertible as a whole), and observe what happens if no row exchange is performed.
(a) For the matrix below, apply the LU decomposition “without row exchanges” directly (set skip_pivot=True),
and observe whether the first elimination step is legitimate.
(b) Decompose the same with the standard PLU (partial pivoting) and compare the two 's. Which one has pivots of larger absolute value? What does this mean for numerical stability?
(c) Find a matrix that is numerically “nearly singular” (determinant close to zero), and observe the last diagonal entry of its .
# Question 4: verification
D = np.array([
[0., 2., 1.],
[3., 4., 0.],
[1., 2., 3.],
], dtype=float)
# (a) Without row exchanges: if D[0,0]=0, the elimination multiplier is inf
print("D[0,0] =", D[0, 0])
print("Without row exchanges, the multiplier of the first elimination step l_{10} = D[1,0]/D[0,0] =",
"∞ (division by zero)" if D[0,0] == 0 else D[1,0]/D[0,0])
steps_D0, _, _, _ = plu_with_steps(D, skip_pivot=True)
print("Last step with skip_pivot=True:", steps_D0[-1]["label"])
print()
# (b) Standard PLU
steps_D, P_D, L_D, U_D = plu_with_steps(D)
print("After the PLU decomposition:")
print("P ="); print(P_D)
print("U ="); print(U_D)
print(f"Diagonal of U (pivots): {np.diag(U_D)}")
print(f"Smallest absolute value of a pivot = {np.min(np.abs(np.diag(U_D))):.4f}")
print()
# (c) Nearly singular matrix
eps = 1e-8
E = np.array([
[1., 2., 3.],
[2., 4., 6.+eps],
[0., 1., 2.],
])
_, _, _, U_E = plu_with_steps(E)
print(f"Last pivot of U for the nearly singular matrix = {U_E[-1,-1]:.2e}")
print(f"Estimated determinant ≈ {np.prod(np.diag(U_E)):.2e}")
Your observations (fill in):
(a)
(b)
(c)
Question 5 (Extension): Determinants and the PLU Decomposition¶
The PLU decomposition provides an efficient way to compute determinants:
where is an matrix, is the number of row exchanges, and the are all diagonal entries (pivots) of .
(a) For Experiment A, compute the number of row exchanges and , and compare with np.linalg.det.
(b) Explain why the lower triangular matrix (with all diagonal entries 1) has determinant 1, so that .
(c) For the rank-deficient matrix (Experiment B), what is ? What does this show?
# Question 5: verification
print("=" * 50)
print(" Determinants and the PLU decomposition")
print("=" * 50)
# Count the row exchanges in Experiment A
swap_count = 0
for s in steps_A:
if "row exchange" in s["label"]:
swap_count += 1
print(f"\nExperiment A: number of row exchanges s = {swap_count}")
pivots = np.diag(U)
prod_pivots = np.prod(pivots)
sign = (-1) ** swap_count
det_plu = sign * prod_pivots
det_numpy = np.linalg.det(A_demo)
print(f"Product of the pivots ∏u_kk = {prod_pivots:.6f}")
print(f"(-1)^s × ∏u_kk = {det_plu:.6f}")
print(f"np.linalg.det(A) = {det_numpy:.6f}")
print(f"The two agree: {np.isclose(det_plu, det_numpy)}")
print()
# Experiment B (rank-deficient)
pivots_B = np.diag(U_B)
print(f"Experiment B, pivots of U: {pivots_B}")
print(f"∏u_kk = {np.prod(pivots_B):.6f}")
print(f"np.linalg.det(B) = {np.linalg.det(B_demo):.6f}")
Your observations (fill in):
(a)
(b)
(c)
Summary of the Experiment¶
After completing the observations above, try to fill in the following blanks in your own words:
| Concept | Your understanding |
|---|---|
| The role of | |
| Why the diagonal entries of are all 1 | |
| The numerical meaning of partial pivoting | |
| Number of rank-one layers = rank of the matrix | |
| The relation between and the pivots |
Hint: the full theoretical explanation of the questions above is in §6.3 of the textbook.