Array Math and solve()

Array Math and solve()

This page picks up where the Arrays page (Arrays) leaves off. There, an array is a list you fill, grow and print; here it is something you compute with as a whole — scale every element in one statement, multiply matrices, and solve systems of equations exactly. Read the Arrays page first; everything below assumes dim, fill, expandable arrays and print of an array.

1. Whole-Array Arithmetic and Reshaping

An expandable array can be assigned a whole array expression: every element is computed at once, and the array takes the shape of the arrays on the right.

dim prices(4)
fill prices with 10, 20, 30, 40
dim with_tax(*)
with_tax = prices * 1.12          ! 11.2, 22.4, 33.6, 44.8
dim discount(4)
fill discount with 1, 2, 3, 4
with_tax = prices - discount      ! array with array: 9, 18, 27, 36
with_tax = with_tax + 5           ! in place: 14, 23, 32, 41
print with_tax                    ! 14, 23, 32, 41

The target must be an expandable REAL array (dim x(*)); it is reshaped to fit whatever the right side produces, so you never size it yourself.

On the right side, any mix of arrays, numbers and + - * / ^ works, including unary minus and parentheses. Every array in one expression must have the same shape, or the catchable exception ARRAYSHAPE is raised. Scalar parts such as rate * 2 are computed once, not once per element.

Shape follows the data. The target of an array assignment takes the shape of the result: a two-dimensional operand makes x two-dimensional, and a later one-dimensional result makes it a plain expandable array again. redim reshapes any array the same way — new bounds or a new number of dimensions — keeping the data in storage order (row by row), and redim x(*) makes an array expandable:

dim x(*)
fill x with 1, 2, 3, 4, 5, 6
redim x(2, 3)                     ! x(1,1)=1 ... x(2,3)=6
redim x(*)                        ! back to one dimension, all six kept
redim x(8)                        ! eight live elements (7 and 8 are zero), still expandable

Because the shape can change while the program runs, a reference with the wrong number of subscripts on a reshaped array (say x(1) while x is two-dimensional) is caught at runtime as the exception WRONGNUMDIMS. An array that is never reshaped keeps its shape for good: for it the subscript count is checked when the program compiles, before anything runs.

Matrix multiply. Two asterisks, **, multiply matrices: with a 2 by 3 and b 3 by 2, z = a ** b makes z 2 by 2 (rows from a, columns from b, each element an exact sum of products). It has the precedence of * and mixes with the other operators: z = a ** b * 2. Both operands must be two-dimensional and the inner sizes must agree, or the catchable exception ARRAYSHAPE names the two shapes. ** is only matrix multiply — between two numbers it is a compile error; exponentiation is ^. A plain vector may stand in for a matrix: on the right of ** it acts as a column, so a matrix times a vector gives a vector; on the left it acts as a row.

dim a(2, 3)
fill a with 1, 2, 3, 4, 5, 6
dim b(3, 2)
fill b with 7, 8, 9, 10, 11, 12
dim z(*)
z = a ** b
print z                           ! 58, 64
                                  ! 139, 154

And you need not assign first just to look: print takes a whole-array expression directly — print z + 1, print v * 2, print a ** b — and prints the result exactly as it prints an array. The same goes wherever a function expects an array: stats$sum(z + 1), stats$mean(v * 2), size(v + 1) all work directly. And the statistics functions do not care about shape — stats$sum(x) on a 2 by 3 array adds all six cells.

2. Solving Equations: solve() and transpose()

A surprising number of everyday questions have the shape "I know the totals — what were the parts?" Each total is one equation; the parts are the unknowns. When there are as many independent totals as unknowns, solve() finds the parts in one statement. Here are three such problems, and what each one teaches.

Problem 1: what does the cafe charge?

A cafe sells only coffee and muffins, and the till records totals, not prices. On Monday it sold 30 coffees and 20 muffins for 130.00; on Tuesday, 25 coffees and 30 muffins for 145.00. What are the two prices?

Why solve(): two unknowns (the prices), two totals (the days). Written out, Monday is 30c + 20m = 130 and Tuesday is 25c + 30m = 145. Put the coefficients in a matrix (one row per day, one column per unknown) and the totals in a vector, and solve(sales, takings) returns the unknowns in column order:

dim sales(2, 2)
fill sales with 30, 20,           ! Monday:  30 coffees, 20 muffins
                25, 30            ! Tuesday: 25 coffees, 30 muffins
dim takings(2)
fill takings with 130, 145        ! the totals: Monday 130.00, Tuesday 145.00
dim price(*)
price = solve(sales, takings)
print price                       ! 2.5, 2.75

Coffee is 2.50 and a muffin 2.75 — printed exactly, not as 2.4999999999999999. solve() uses fraction-free elimination in Sheerpower's exact decimal arithmetic: the only place an answer can round is one final division per unknown, so an answer that terminates comes out exact.

Problem 2: how many of each vehicle?

A delivery company knows it runs 8 vehicles, that together they carry 21 tons, and that the fleet costs 1650 a day. A van carries 1 ton and costs 100 a day; a truck 2 tons and 150; a lorry 5 tons and 400. How many vans, trucks and lorries are there?

Why solve(): three unknowns, three independent facts — a count, a capacity, a cost. Each fact is a row: what one van, one truck and one lorry contribute to it.

dim fleet(3, 3)
fill fleet with 1,   1,   1,      ! vehicles:  one each
                1,   2,   5,      ! tons:      van 1, truck 2, lorry 5
              100, 150, 400       ! cost/day:  van 100, truck 150, lorry 400
dim facts(3)
fill facts with 8, 21, 1650      ! the totals: 8 vehicles, 21 tons, 1650 a day
dim count(*)
count = solve(fleet, facts)
print count                       ! 1, 5, 2

One van, five trucks, two lorries. The same shape serves any "blend" question — ingredients meeting nutritional targets, investments meeting a return and a risk figure, staffing meeting hours and budget.

Problem 3: which way are sales heading?

Six months of sales: 13, 13, 15, 19, 20, 22. They rise, but not in a straight line — no line passes through all six points. The question is the best straight line, the one that misses the points by the least (in the least-squares sense), and what it predicts for month 7.

Why transpose(): this time there are more facts (six months) than unknowns (a slope and an intercept), so the system cannot be solved as it stands. The classic remedy is to multiply both sides by the transpose of the coefficient matrix, which produces a square system — the normal equations — whose solution is the best-fit line. Build the coefficient matrix with one row per month (month, 1), and the two lines of algebra become one line of Sheerpower. You do not need the theory to use it: the pattern solve(transpose(a) ** a, transpose(a) ** y) is the whole trick, and it fits a line (or a plane, with more columns) through any set of points.

dim design(6, 2)                  ! one row per month: (month, 1)
fill design with 1, 1,  2, 1,  3, 1,  4, 1,  5, 1,  6, 1
                                  ! month multiplies the slope, the 1 the intercept
dim sales(6)
fill sales with 13, 13, 15, 19, 20, 22   ! the six monthly figures, same order
dim fit(*)
fit = solve(transpose(design) ** design, transpose(design) ** sales)
print fit                         ! 2, 10
print fit(1) * 7 + fit(2)         ! 24: fit(1) is the slope, fit(2) the intercept

Sales grow by 2 a month from a base of 10, and month 7 should bring 24. (transpose(design) ** sales multiplies a matrix by the plain sales vector — the column rule from section 1.)

One honest expectation to set: this data was chosen so the line comes out to exactly 2 and 10. Real best-fit lines almost never terminate — a slope of 1.8571428571428571 is the normal case — and a long decimal is the arithmetic being honest, not wrong. The rule from Problem 1 still holds: an answer that terminates comes out exact, one that does not is rounded once, at the end.

transpose() on its own: rows into columns

Data often arrives one way and is wanted the other. A shop sells coffee beans and tea, and its sales come in as one row per month with a column per product: in month 1 it sold 120 bags of beans and 80 of tea, in month 2 150 and 95, in month 3 170 and 110. A report wants one row per product with the months across. transpose() turns the 3 by 2 into a 2 by 3:

dim by_month(3, 2)                ! one row per month: (beans, tea)
fill by_month with 120, 80,       ! month 1
                   150, 95,       ! month 2
                   170, 110       ! month 3
dim by_product(*)
by_product = transpose(by_month)
print by_product                  ! 120, 150, 170   -- beans, months 1 to 3
                                  ! 80, 95, 110     -- tea

When there is no answer

If the facts are not independent — Tuesday was simply double Monday — the system has no unique solution and solve() raises the catchable exception SINGULAR. Catch it where a bad data set is a real possibility (exception handling has a page of its own: Exception Handling):

when exception in
  price = solve(sales, takings)
use
  if extype = exceptiontype('singular') then
    print 'these totals do not pin down the prices'
  end if
end when

Shape mistakes — a matrix that is not square, a totals vector with the wrong number of rows — raise ARRAYSHAPE, and the message names the shapes it saw.

(Show/Hide Array Math Takeaways)