Start with linalg: and write one statement per line: matrices as A = [1 2; 3 4], column vectors as v = (1, 2), combined with +, -, *, ^T and ^-1.
Functions such as det(A), rank(A) or eigenvalues(A) return results that can be assigned and reused, e.g. d = det(A) + rank(A). Capitalize the first letter (Det(A)) to also get a step-by-step explanation of the calculation.
Annotate an operator with a unique marker to animate it: C = A * @product B. In the declarations, operation=product; step=0; selects the marked operation and the step to show (step=-1 is the final result). The Animate ... conversion buttons turn the calculation into a click-through animation.
linalg:
A = [1 2; 3 4]
B = [0 1; 1 0]
d = Det(A)
C = A * @product B
--declarations--
operation=product;
step=0;
--declarations-end--
Animations, declarations and further functions
Optionally mark an operator: C = A * @product B. Select operation=product; step=0; in declarations to show its initially empty result. Each step adds one scalar product to the running sum (m*n*k steps); unfilled cells remain distinct from zeros. Use step=-1 for the final result or the Animate multiplication conversion for a click-through SVG. Marker names must be unique; quote names with spaces, e.g. @"my product", and URL-encode them in declarations (operation=my+product;). For MLP training animations, define input x, target, weights and biases, then declare mlpLayers=#W1,b1,tanh;W2,b2,sigmoid#; mlpStep=0;. MLP targets accept a scalar (target=0.8) or row/column vector (e.g. target=[0.3,0.7]); they are normalized to columns internally. Animate MLP training shows a network diagram and the current calculation's matrices, one neuron per step, through forward propagation, backpropagation, one SGD update and another forward pass comparing old and new activations. Each frame explains its calculation; the loss frame compares prediction and target and explains mean squared error or cross-entropy. For a transformer self-attention animation, define token rows X and weights Wq, Wk, Wv (optionally Wo), then declare attention=#X,Wq,Wk,Wv#; attnTokens=#the,cat,sat#; attnStep=0; (attnCausal=true; masks later tokens). Animate attention steps through queries, keys, values, scaled scores, row-wise softmax weights and the weighted sum of values, with a token graph whose line thickness shows the attention weights. Conversion buttons are shown only when applicable: MLP controls need MLP declarations, multiplication controls need a marked matrix product, and Next is hidden at the final step. Use svgHeight=800; in declarations to start with a taller web output box (pixels; 0 keeps the default). This works for all diagram types; manual resizing remains available. The buttons above the drawing zoom out/in, fit it to the width, or fit the whole drawing into the box. Standalone calls such as det(A), rank(A), span(A), eigenvalues(A) and eigenvectors(A) return reusable results. DesignMatrix(x, degree) builds the polynomial-regression matrix $(1,x_i,\ldots,x_i^{degree})$ (degree 0 through 12); use showPlot=true to plot each power feature curve with the sample points. DesignMatrix(x, degree, y) also fits the least-squares polynomial to the responses y and plots it (blue) with the data points (red) and residuals. svd(A) and evd(A) display their decompositions in place. Use mean(data), variance(data) and std(data) for statistics; capitalize them for worked steps. Classification and probability use precision(tp, fp), recall(tp, fn) and bayes(prior, likelihood, evidence). binomial(n, p, k) and poisson(lambda, k) return P(X=k). Neural networks use dense(W, X, b), relu(Z), sigmoid(Z), tanh(Z), softmax(Z) and explicit backward functions. Capitalized calls teach the concepts step by step; showPlot=true also draws activation/slope curves and softmax class shares. Set showPlot=true in the declarations to draw 2D/3D objects and empirical histograms with mean and standard-deviation markers.
Neural networks and explicit backpropagation
Columns are examples; rows are features, neurons or classes. dense(W, X, b) computes Z=W*X+b with shapes W[o,i], X[i,B], b[o,1], Z[o,B]. Only this operation broadcasts the bias column. For one example use a column vector; [x] is a one-feature example. A softmax row vector means one class and multiple examples. Ordinary addition still requires equal shapes.
relu(Z), sigmoid(Z) and tanh(Z) operate elementwise on scalars or matrices. softmax(Z) normalizes each matrix column stably. hadamard(A, B) multiplies matching entries, unlike matrix multiplication. mse(A, T) averages squared errors over all entries, without a factor 1/2. crossentropy(Z, T) takes logits, not probabilities: it fuses softmax and cross-entropy with stable log-sum-exp and averages over example columns. Target columns must be probabilities summing to one (tolerance 1e-12); one-hot and soft labels are supported. The actual target sum is retained when differentiating.
Start backpropagation with msegrad(A, T) or crossentropygrad(Z, T). These already include the loss averaging factor. reluback(Z, GA), sigmoidback(Z, GA), tanhback(Z, GA) and softmaxback(Z, GA) take the original preactivation/logits Z and upstream loss gradient GA, returning GZ. ReLU uses derivative zero at zero; softmaxback includes cross-class coupling. Do not apply softmaxback after crossentropygrad: the fused gradient already includes it.
densegradw(GZ, X) returns GZ*X^T; densegradb(GZ) sums example columns; densegradx(W, GZ) returns W^T*GZ for the preceding layer. No extra averaging is performed. Save forward inputs and preactivations, compute all gradients using original weights, then update Wnew=W-eta*gW and bnew=b-eta*gb. Operations return ordinary reusable values, not hidden network objects; there is no automatic differentiation tape.
linalg:
W=[0.2 0.3]
x=(1,2)
b=[0.1]
target=[1]
z=Dense(W,x,b)
a=Sigmoid(z)
loss=Mse(a,target)
ga=Msegrad(a,target)
gz=Sigmoidback(z,ga)
gW=Densegradw(gz,x)
gb=Densegradb(gz)
gx=Densegradx(W,gz)
Wnew=W-0.1*gW
bnew=b-0.1*gb
Capitalized variants show real per-entry calculations, derivative conventions and ML uses. Named examples include a two-layer network and batched classification with loss-reducing SGD updates. Nonfinite entries, empty inputs, shape mismatches and overflowing intermediates are explicit errors. Sigmoid/tanh saturation and tiny gradients may still underflow in double precision. Neural operations do not add architecture plots.
Beginner explanations introduce weighted votes and bias, preactivation versus activation, targets and losses, local slopes, upstream loss gradients and the chain rule. Backpropagation computes sensitivities; gradient descent separately updates parameters with a learning rate. Large steps need not decrease loss.
With showPlot=true, ReLU, sigmoid and tanh calls show activation and local-slope curves side by side over [-5,5], with separate vertical scales and red markers for up to 12 distinct supplied preactivations. Off-window and capped marker counts are reported. ReLU's derivative at zero uses a filled dot at the chosen zero and an open dot at one. Slopes are not weight updates or upstream-multiplied gradients. Softmax is coupled: it shows probability bars, not a scalar activation curve. At most four example columns and 12 classes per column are shown, with omitted class mass and limits disclosed; no renormalization. Fused cross-entropy contributes its logits' softmax view. Standalone, nested, lowercase and backward calls also plot; repeated views share plots and input snapshots survive redefinitions. No training history is inferred. The named ANN-basics example demonstrates these views. Use showFormulas=false; for plots only.
Future convolution and pooling can follow the same explicit forward/backward contract. Convolution, pooling, multichannel tensors, automatic differentiation and training loops are not implemented yet; channel layout will need an explicit definition.
Binomial and Poisson distributions
binomial(n, p, k) returns P(X=k) for n independent trials with constant success probability p: choose(n,k)*p^k*(1-p)^(n-k). poisson(lambda, k) returns exp(-lambda)*lambda^k/k!: an event-count model with expected count lambda in a fixed interval. A constant event rate times exposure gives lambda. These return scalar PMFs, not cumulative probabilities or samples. Use Binomial(10, 0.3, 3) or Poisson(4, 2) for substituted formulas, stable log-probability steps, mean, variance and ML uses.
All inputs must be finite scalars. n and k must be nonnegative integers at most 9007199254740991 (the exact integer range used here); p is in [0,1], and lambda is in [0,9007199254740991]. Binomial k>n returns exact zero; negative or fractional counts are errors. n=0, p=0/1 and lambda=0 are handled explicitly. A nonzero requested probability that underflows double precision is an explicit error.
With showPlot=true, each distribution call adds a discrete PMF plot, also inside assignments or nested expressions. Repeated identical calls share one plot; red highlights k, or the caption states that k lies outside the window. Binomial n<=100 shows full support; larger binomials and Poisson show up to 101 consecutive outcomes centered near the mean (approximately six standard deviations plus six outcomes on either side, capped at fifty). Captions report the truncated window and estimated probability outside it. Masses are not renormalized; underflowed plotted tail masses are explicitly reported as drawn zeros.
Precision, recall and Bayes' formula
precision(tp, fp) computes TP/(TP+FP): the fraction of predicted positives that are correct. recall(tp, fn) computes TP/(TP+FN): the fraction of actual positives found (sensitivity). TP, FP and FN mean true positives, false positives and false negatives. Arguments must be finite nonnegative scalar counts; fractional counts are allowed for weighted observations. A zero denominator is undefined and raises a line-numbered error, not zero or one.
bayes(prior, likelihood, evidence) computes P(H|E)=P(H)*P(E|H)/P(E). Supply probabilities in [0,1], with evidence P(E)>0; evidence is not a false-positive rate. Consistent inputs require prior*likelihood ≤ evidence ≤ prior*likelihood+(1-prior). For a diagnostic test with prior=0.01, sensitivity=0.9 and false-positive rate=0.05, first compute evidence=0.9*0.01+0.05*(1-0.01), then Bayes(0.01, 0.9, evidence), giving about 0.1538, not 0.9.
Precision(tp, fp), Recall(tp, fn) and Bayes(prior, likelihood, evidence) show definitions, substituted inputs, denominator or joint-probability calculation, final ratio and ML uses. All return reusable scalars and can be nested in matrix expressions. They do not produce histograms: counts and probabilities are not observation vectors. High-precision decimal intermediates avoid premature overflow/underflow; unrepresentable nonzero results are explicit errors.
Statistics in the same script
Mix statistics and matrix algebra freely: data=[2 4 4 4 5 5 7 9], mu=mean(data), v=variance(data), s=std(data) return scalars 5, 4 and 2. s*A scales a matrix; mean(eigenvalues(A)) composes functions. Mean(data), Variance(data) and Std(data) show running sums, deviations, squared deviations and ML uses as appropriate.
Each finite entry of a nonempty row or column vector is one equally weighted observation. Use [2] for a singleton; scalars, empty vectors and general matrices are rejected instead of guessing an aggregation axis. Variance is population variance, sum of squared deviations divided by n, not n-1. A singleton has zero variance and std. High-precision decimal intermediates avoid premature overflow; an unrepresentable requested scalar result is an explicit error.
With showPlot=true, statistical calls also produce empirical relative-frequency histograms, including nested or assigned calls. Identical data share one histogram. The bin count is ceil(sqrt(n)), capped at 12; constant data use one bin. Red marks the mean and dashed lines mark mean plus/minus one standard deviation. Bin masses (not bar areas) sum to one. No normal distribution or fixed coverage probability is assumed; the shaded interval is not a confidence interval.
Syntax and expressions
Names are case-sensitive, start with a letter and may contain letters, digits and underscores. Define names before using them. Spaces (or commas) separate matrix entries; semicolons separate rows. A vector (1, 2, 3) is a column vector; [1 2 3] is a row vector. Decimal points and scientific notation (1e-3) are supported. // starts a line comment.
Expressions support +, -, *, unary signs, parentheses, ^T (transpose) and ^-1 (inverse). Multiplication accepts scalars and dimension-compatible matrices or vectors. Addition and subtraction require equal dimensions. Expressions remain in the script; only the output shows their evaluated results.
Functions
det(A): determinant of a square matrix.rank(A): number of independent columns.trace(A): diagonal sum of a square matrix.norm(A): Euclidean vector norm or Frobenius matrix norm.rref(A): reduced row echelon form.basis(A): independent original columns as a basis matrix.nullspace(A): kernel basis, with one basis vector per column. A trivial kernel has an n-by-0 basis matrix.transpose(A) and inv(A): alternatives to A^T and A^-1. Inversion requires a nonsingular square matrix.
These functions work on their own line or inside expressions, for example d = det(A) + rank(A) or R = rref(A). Rank and basis calculations use a relative numerical tolerance of 1e-12.
Lowercase det calls of size two or larger display one first-row Laplace expansion: det(A) = a11 * det(M11) - a12 * det(M12) + ... = result. Each minor is shown as a matrix with the first row and the corresponding column removed; its determinant is not expanded further. A 1x1 determinant shows its entry directly. Traces up to size four and norms with up to four entries also show the substituted formula. This also works for d = det(A); nested calls remain concise.
eigenvalues(A) returns a column of real eigenvalues in descending order, including algebraic multiplicities. eigenvectors(A) returns unit-length eigenvectors as columns in that eigenvalue order. Both work in assignments and expressions. For diagonalizable matrices the entries and columns correspond one-to-one. For defective matrices only independent eigenvectors per eigenspace are returned, so there are fewer columns and no entry-by-entry pairing with the full eigenvalue list. Signs and bases for repeated eigenvalues are not unique; nonsymmetric eigenvectors need not be orthogonal. Both require a square matrix and entirely real spectrum; invalid inputs and numerical failures report errors. Eigenvalues and Eigenvectors explain their meaning and algorithm. eig/Eig remain compatibility aliases for evd/Evd.
A = [4 1; 2 3]
values = eigenvalues(A)
vectors = eigenvectors(A)
Evd(A)
span(A) displays the span of A's columns using independent generators. span(u, v, ...) accepts column vectors of equal dimension. svd(A) displays the singular value decomposition; evd(A) displays a real eigenvalue decomposition, or explains why it is unavailable. These three calls must stand on their own line and cannot be assigned to a name. Use basis(A) for a reusable matrix representing the column space.
Capitalize the first letter for an algorithm explanation before the result: Det(A), Rank(A), Rref(A), Basis(A), Nullspace(A), Inv(A), Trace(A), Norm(A), Transpose(A), Span(A), Svd(A), Evd(A), Eigenvalues(A), Eigenvectors(A). Elimination-based functions show row operations and intermediate matrices. Nullspace(A) additionally explains Ax=0, rank-nullity, free variables, pivot equations, each basis-vector construction and the general solution x=Nt. It distinguishes the trivial kernel from its empty basis. Svd(A) shows actual Jacobi rotations, working matrices after each sweep, column lengths and sorted normalized factors. Both SVD forms use the thin decomposition: for m rows, n columns and k=min(m,n), U is m-by-k, Sigma is k-by-k and V is n-by-k. The displayed V* is V-transpose (k-by-n). A full SVD has square U and V with extra directions multiplying zeros in a rectangular Sigma. Svd explains that V-transpose swaps rows and columns, while V* is the conjugate transpose (also V^H), additionally replacing a+bi by a-bi. They are equal for the real matrices supported here. V* takes inner products with the right singular vectors; V*V=I, but rectangular thin V has no two-sided inverse. Det(A) omits the Laplace expansion and shows the elimination followed by the result. Other capitalized functions also show input-specific work: running diagonal sums for Trace, stable hypot accumulation for Norm, row-to-column moves for Transpose, and original-column selection for Basis/Span. Evd, Eigenvalues and Eigenvectors show Jacobi or QR evolution, spectral blocks and eigenspace construction. Long runs show the first eight Jacobi sweeps and every tenth, or first eight QR iterations and every hundredth, plus the final spectral matrix when reached; selected snapshots are labeled. Differential operators show differentiation rules, intermediate partials, component assembly and point substitution. Every capitalized function, including the differential operators below, also summarizes its use in machine learning (ML). Results and source stay unchanged; numeric calls also work in expressions such as d = Det(A) + Rank(A).
Example
linalg:
A = [1 2 3; 2 4 6; 1 1 1]
rank(A)
span(A)
N = nullspace(A)
Z = A * N
svd(A)
Functions and differential operators
f(x, y) = x^2 + 2*x*y + 3*y^2
p = (1, 2)
g = grad(f, p)
Grad(f, p)
grad(f, p) returns the gradient as a column vector (here, (6, 14)). Grad(f, p) first displays the symbolic partial derivatives in parameter order and the evaluation point. Derivatives are symbolic, not finite differences; gradient calls also work in matrix expressions.
With showPlot=true, direct grad/Grad results show a local function curve with point and tangent for one variable, or a contour map with labeled levels and a unit gradient-direction arrow for two variables. The window extends two units on either side of the point. Undefined samples are omitted with a notice; sampled plots may miss discontinuities. Aliases retain the view, transformed vectors do not. Three-variable functions keep the existing 3D gradient arrows without level surfaces. Original vector plots retain true magnitudes.
Scalar-function bodies support + - * / ^, parentheses, sin, cos, exp, and natural logarithm log. Vector-valued functions use tuples, e.g. F(x, y) = (x*y, sin(x)). Use only declared parameters and numeric constants; references to scene objects and calls to other user-defined functions are not supported. Define functions before use. Powers associate to the right; -x^2 means -(x^2). Direct function evaluation is not yet supported.
jacobian(F, p): m-by-n matrix of first derivatives; rows are output components, columns are input parameters. Scalar functions have a one-row Jacobian.hessian(f, p): n-by-n matrix of second derivatives of a scalar function in parameter order; symmetric when the second derivatives are continuous.laplacian(f, p): sum of pure second derivatives, the trace of the Hessian.divergence(F, p): sum of derivatives of component i with respect to variable i; requires one output per input. This is the trace of the Jacobian.curl(F, p): three-component curl; requires three inputs and three outputs.
All operators use symbolic derivatives at a column-vector point. Jacobian, Hessian, Laplacian, Divergence and Curl explain the operator and show symbolic derivatives before the result. Results are reusable in expressions. Curl arrows start at the evaluation point. Jacobians and Hessians use matrix plots.
f(x, y) = x^2 + 2*x*y + 3*y^2
p = (1, 2)
Hessian(f, p)
Laplacian(f, p)
F(x, y) = (x*y, sin(x))
Jacobian(F, p)
Divergence(F, p)
G(x, y, z) = (-y, x, 0)
Curl(G, (1, 2, 3))
The evaluation point must be a column vector with one entry per parameter; for a one-variable function, a scalar such as p = (2) or 2 is also accepted, as is [2]. Parentheses still group scalar arithmetic; the operator interprets a scalar as a one-coordinate point. Invalid domains and non-finite derivatives produce line-numbered errors. With plotting enabled, assigned gradient vectors start at their evaluation point, not at the origin. Matrix operations on them produce ordinary vectors drawn from the origin. Standalone gradient calls are not plotted; contours and gradient fields are not yet available.
Plots and declarations
The following block enables plots. Set showFormulas=false to show only the plot. Formulas are always shown when plotting is disabled.
--declarations--
showPlot=true;
showFormulas=true;
--declarations-end--
2D/3D column vectors are drawn as arrows. A 2x2 matrix maps the unit square; a 3x3 matrix maps the unit cube, shown with shaded faces in a fixed oblique view. Mixed scenes produce separate 2D and 3D pictures. Other dimensions still support calculations but are not plotted. For a basis change, put the new basis vectors in the columns of B and write c = B^-1 * v.
Legacy matrix: scripts with spaces between columns and newlines between rows are still accepted. Invalid inputs report their line number.