Matrix Inversion by DeepSeek

Discussions related to mathematics, numerical methods, graph plotting etc.
Richard Russell
Posts: 272
Joined: Tue 18 Jun 2024, 09:32

Matrix Inversion by DeepSeek

Post by Richard Russell »

I've never been sure whether the Matrix Inversion procedure in arraylib.bbc is any good, so I thought i'd ask DeepSeek to write a new one so I could compare them. Here is the question I asked:
Please write a BBC BASIC procedure to invert an arbitrary-size square matrix. The input matrix will be in a 2-dimensional array and the output matrix may either be in the same array or a different array, remember that BBC BASIC array indexes are zero-based. You may use BBC BASIC 5 features, including array arithmetic, if they are helpful. You should prioritize speed over code size.
Based on our recent previous experience I expected DeepSeek to make working code first time, but it didn't: there were syntax errors resulting from it not understanding BBC BASIC properly (most obviously that it seemed to think BBC BASIC has array-slicing functionality; I wish that it did but...).

It was easy enough to correct its mistakes but it's a bit of a black mark that it suffered from some hallucination on this task. Anyway with the corrections it runs and gives near enough the same answers and takes near enough the same time as arraylib, so I'm not able to say which is 'better'.

But as always what I found most fascinating was not the end result but its 'thinking'. Is there genuine reasoning here or not, and where did it get that stuff about slicing and assigning entire rows?:
Okay, I need to write a BBC BASIC procedure to invert a square matrix of arbitrary size. Let me think about how to approach this. The input is a 2D array, and the output can be in the same array or a different one. The indexes are zero-based. Since BBC BASIC 5 allows array arithmetic, maybe that can help, but I should prioritize speed over code size.

First, matrix inversion is typically done using Gaussian elimination, often with partial pivoting to avoid division by zero and reduce numerical errors. The steps are: form an augmented matrix with the identity matrix, perform row operations to convert the original matrix to the identity, and the augmented part becomes the inverse. So I need to implement that.

Let me outline the steps:

1. Check if the matrix is invertible. If during the process a pivot is zero even after row swaps, it's singular. So I need to handle that, maybe by returning an error?

But the problem statement doesn't mention error handling, just to write the procedure. So perhaps assume the input is invertible.

2. Create an augmented matrix. The original matrix and the identity matrix. Then perform row operations.

But in BBC BASIC, how are 2D arrays handled? Let's say the input is DIM A(n%,n%), then we need to create a bigger matrix for the augmented part. Wait, the input is a square matrix of size N x N. So the augmented matrix is N x 2N.

But in code, perhaps create a working matrix that's N x 2N. Then after performing the row operations, the right half is the inverse. Then copy that back into the output array.

Wait, but the user may want the inverse to be in the same array or a different one. So perhaps the procedure should take two parameters: the input array and the output array. The output array can be the same as input, but perhaps we need to create a temporary augmented matrix.

Alternatively, since BBC BASIC allows array operations, maybe we can perform in-place inversion. But I think the standard method requires the augmented matrix. So here's a plan:

- The procedure will take two arrays: the input matrix and the output matrix. They can be the same.

- The steps:

a) Get the size of the matrix. Since it's square, check that the input array is square.

But the code may need to get N% from the input array. In BBC BASIC, the DIM function returns the size, but since arrays can have different lower bounds, perhaps using DIM(input(),1) to get the second dimension's upper bound. Wait, BBC BASIC arrays are zero-based here. Wait, the user says the indexes are zero-based. So if the matrix is N x N, then it's DIM A(N-1, N-1). Wait, no. If it's zero-based, then for an N x N matrix, the maximum index is N-1 in each dimension. So the size is (N-1)+1 = N. So for example, a 3x3 matrix would have indexes 0 to 2 in each dimension.

So to get the size, perhaps the code can get the second dimension's size. So if the input is a two-dimensional array, then N% = DIM(input(), 1) + 1. Wait, in BBC BASIC, DIM(array(), dimension) returns the upper bound. For a zero-based array, the upper bound for each dimension is N-1. So for example, a 3x3 matrix is DIM A(2,2). So DIM(A(), 0) returns 2 (the upper bound of the first dimension), same for DIM(A(),1). So the size N% is DIM(A(),0)+1.

So first, get N% as DIM(input(),0) +1. Then check that the input is a square matrix by ensuring DIM(input(),1) is the same as DIM(input(),0).

Once N% is determined, create the augmented matrix. Then perform Gaussian elimination with partial pivoting.

But how to handle this in code. Let's outline:

- For each column (from 0 to N-1):

- Find the pivot row: the row with the maximum absolute value in the current column, starting from the current row.

- If the pivot is zero, the matrix is singular (but we proceed as per user's assumption?)

- Swap the current row with the pivot row in the augmented matrix.

- Normalize the pivot row so that the pivot element is 1.

- Eliminate all other rows by subtracting the appropriate multiple of the pivot row.

Wait, the standard algorithm for Gaussian elimination for inverse:

1. Augment the matrix with the identity matrix.

2. For each column (which is the current pivot column):

a. Find the row with the maximum element in that column (from current row downwards) to choose the pivot.

b. Swap the current row with the pivot row.

c. Divide the current row by the pivot element to set the pivot to 1.

d. For all other rows, subtract the current row multiplied by the element in their column (of the current pivot column) to eliminate that element.

But in code, this would require nested loops.

So, the steps in code:

Given the input matrix, we need to create an augmented matrix. Let's say the original matrix is A and the augmented part is the identity. So for each row i, columns 0 to N-1 are A's elements, columns N to 2N-1 are the identity.

But in BBC BASIC, how to handle that. The code could create a temporary matrix augmented which is N x 2N.

But perhaps in the code, the input and output are separate arrays. Wait, the user says that the output may be in the same array or a different one. So perhaps the procedure is:

DEF PROCinvert(input(), output())

But perhaps the input and output can be the same array. So during processing, the code would need to create a temporary augmented matrix. So the steps would be:

1. Check if input is square. Get N% as DIM(input(),0)+1. Check that DIM(input(),1) = DIM(input(),0). If not, error.

2. Create a temporary augmented matrix, aug(), of size N% x 2*N%.

3. Initialize the left half of aug() with input(), and the right half with the identity matrix.

4. Perform Gaussian elimination on aug().

5. After elimination, the right half of aug() is the inverse. Copy this into output().

But BBC BASIC allows array operations, perhaps. So perhaps for steps 3 and 5, there are efficient ways to copy.

But the code has to be written in BBC BASIC. Let's think about how to code this.

First, the code outline:

DEF PROCinvert(input(), output())

LOCAL N%, aug(), i%, j%, k%, pivot_row%, max_val, temp, factor

N% = DIM(input(), 0) + 1

REM Check if the input is square

IF DIM(input(), 1) + 1 <> N% THEN ERROR 100, "Matrix is not square"

REM Create augmented matrix: N% rows, 2*N% columns

DIM aug(N%-1, 2*N% -1)

REM Copy input into left half, identity into right half

FOR i% = 0 TO N% -1

FOR j% = 0 TO N% -1

aug(i%, j%) = input(i%, j%)

NEXT

FOR j% = N% TO 2*N% -1

aug(i%, j%) = 0

NEXT

aug(i%, i% + N%) = 1

NEXT

REM Perform Gaussian elimination with partial pivoting

FOR k% = 0 TO N% -1

REM Find pivot row

pivot_row% = k%

max_val = ABS(aug(k%, k%))

FOR i% = k% +1 TO N% -1

IF ABS(aug(i%, k%)) > max_val THEN

pivot_row% = i%

max_val = ABS(aug(i%, k%))

ENDIF

NEXT i%

REM If max_val is zero, matrix is singular (handle error?)

IF max_val = 0 THEN ERROR 101, "Matrix is singular"

REM Swap rows k% and pivot_row%

IF pivot_row% <> k% THEN

FOR j% = 0 TO 2*N% -1

temp = aug(k%, j%)

aug(k%, j%) = aug(pivot_row%, j%)

aug(pivot_row%, j%) = temp

NEXT j%

ENDIF

REM Normalize the pivot row

temp = aug(k%, k%)

FOR j% = k% TO 2*N% -1

aug(k%, j%) = aug(k%, j%) / temp

NEXT j%

REM Eliminate other rows

FOR i% = 0 TO N% -1

IF i% <> k% THEN

factor = aug(i%, k%)

FOR j% = k% TO 2*N% -1

aug(i%, j%) = aug(i%, j%) - factor * aug(k%, j%)

NEXT j%

ENDIF

NEXT i%

NEXT k%

REM Copy the right half of aug() to output()

FOR i% = 0 TO N% -1

FOR j% = 0 TO N% -1

output(i%, j%) = aug(i%, j% + N%)

NEXT j%

NEXT i%

ENDPROC

Wait, but in the elimination step, after normalizing the pivot row, for each other row, subtract factor times the pivot row. But the loop for j% can start from k%, since elements before k% in the row are already zero.

But in code, the inner loop for j% runs from k% to 2*N%-1. That's correct.

But when normalizing the pivot row, we divide all elements from column k% to 2*N%-1. The elements before k% are already zero, except for the identity part. Wait, no. Wait, in the augmented matrix, during the process, the left part is transformed into the identity, and the right part becomes the inverse. But the code should handle the elimination correctly.

Wait, but in the code above, during the elimination step for each row i% != k%, it subtracts factor * pivot row from row i%. But factor is aug(i%,k%) which, after the pivot row has been normalized, is the element to eliminate. So after the subtraction, aug(i%,k%) becomes zero.

Yes.

But after all steps, the left half of aug() becomes the identity matrix, and the right half is the inverse.

But wait, when we process each column k%, the pivot is in row k%, and after normalization, we eliminate all other rows in the column k%. So that's correct.

Now, in BBC BASIC, array arithmetic is allowed. But in the code above, loops are used. However, using array operations could make the code faster. But for arbitrary size matrices, perhaps the loops are necessary. But given that the user says to prioritize speed, perhaps using array slices where possible.

But BBC BASIC's array arithmetic may allow for vector operations. For example, after swapping the rows, instead of looping through each j%, can we swap the entire rows? Like, in BBC BASIC, you can assign a whole row with another. For example, if you have two rows, row1 and row2, you can do SWAP row1, row2. Wait, but in BBC BASIC, arrays are handled as individual elements. So perhaps not. So the code as written is correct.

But let me check.

In BBC BASIC, you can perform array operations. For example, if you have two arrays A() and B(), you can do A() = B() * 2. But when dealing with rows, perhaps you can take a slice.

Wait, BBC BASIC allows you to use a 'slice' of an array. For example, to access a row, you can use aug(i%,) which represents the entire row. But I'm not sure if that's the case. Let me recall.

In BBC BASIC, the comma is used to separate dimensions. So to refer to an entire row, you can leave the second dimension unspecified. For example, A(i%,) refers to the i%th row as a vector. So you can perform operations like:

A(i%,) = A(i%,) / temp

Which would divide all elements in row i% by temp. Similarly, you can assign a row to another row.

So perhaps, in the code, instead of looping through j%, we can use these vector operations for speed.

Similarly, when swapping rows, we can swap the entire rows.

So let's revise the code to take advantage of this.

For example:

When swapping rows:

Instead of looping through all j% to swap each element, perhaps:

temp_row() = aug(k%,)

aug(k%,) = aug(pivot_row%,)

aug(pivot_row%,) = temp_row()

But wait, in BBC BASIC, can you assign a whole row like that?

Yes, in BBC BASIC, you can assign an entire row to another row if the dimensions match. So perhaps using this would be faster than looping through each element.

Similarly, when normalizing the pivot row:

temp = aug(k%,k%)

aug(k%,) = aug(k%,) / temp

But wait, in BBC BASIC, when you write aug(k%,) = aug(k%,) / temp, this is an array operation. Each element in row k% is divided by temp. So that's much more efficient than looping through j%.

Similarly, during elimination:

factor = aug(i%,k%)

aug(i%,) = aug(i%,) - factor * aug(k%,)

This would subtract factor times the pivot row from row i%. But note that the columns here start from k%, but the code above would do it for all columns. But since in the earlier steps, the pivot row has been normalized and the previous columns have been processed, the elements before column k% are already zero in the pivot column (except in the pivot row), so subtracting in those columns would have no effect, but it's okay to do so. However, in the code, when we normalize the pivot row, we only start from column k%. Wait, in the code above, during normalization, the loop runs from j% = k% to 2*N% -1. So if we can use array operations starting from column k%, then perhaps we can optimize.

But in BBC BASIC, array slices can have ranges. For example, aug(k%,k% TO 2*N%-1) would represent the elements from column k% to 2*N%-1 in row k%.

So maybe:

When normalizing, instead of looping j% from k% to 2*N% -1:

aug(k%,k% TO 2*N%-1) = aug(k%,k% TO 2*N%-1) / temp

Similarly, during elimination:

aug(i%,k% TO 2*N%-1) = aug(i%,k% TO 2*N%-1) - factor * aug(k%,k% TO 2*N%-1)

That would be more efficient, as it avoids looping through each element in BASIC code and uses built-in array operations which are faster.

So revising the code:

DEF PROCinvert(input(), output())

LOCAL N%, aug(), i%, j%, k%, pivot_row%, max_val, temp, factor, temp_row()

N% = DIM(input(), 0) + 1

REM Check if the input is square

IF DIM(input(), 1) + 1 <> N% THEN ERROR 100, "Matrix is not square"

DIM aug(N%-1, 2*N% -1)

REM Copy input into left half and identity into right half

FOR i% = 0 TO N% -1

aug(i%, 0 TO N%-1) = input(i%,) COPY input row into left half

aug(i%, N% TO 2*N%-1) = 0 Initialize right half to zero

aug(i%, i% + N%) = 1 Set identity element

NEXT

FOR k% = 0 TO N% -1

REM Find pivot row

pivot_row% = k%

max_val = ABS(aug(k%, k%))

FOR i% = k% +1 TO N% -1

current_val = ABS(aug(i%, k%))

IF current_val > max_val THEN

pivot_row% = i%

max_val = current_val

ENDIF

NEXT i%

IF max_val = 0 THEN ERROR 101, "Matrix is singular"

REM Swap rows k% and pivot_row%

IF pivot_row% <> k% THEN

SWAP aug(k%,), aug(pivot_row%,)

ENDIF

REM Normalize the pivot row from column k% onwards

temp = aug(k%, k%)

aug(k%, k% TO 2*N%-1) = aug(k%, k% TO 2*N%-1) / temp

REM Eliminate other rows

FOR i% = 0 TO N% -1

IF i% <> k% THEN

factor = aug(i%, k%)

aug(i%, k% TO 2*N%-1) = aug(i%, k% TO 2*N%-1) - factor * aug(k%, k% TO 2*N%-1)

ENDIF

NEXT i%

NEXT k%

REM Copy the inverse from the right half of aug() to output()

FOR i% = 0 TO N% -1

output(i%,) = aug(i%, N% TO 2*N%-1)

NEXT i%

ENDPROC

Wait, but in BBC BASIC, can we assign a slice like N% TO 2*N% -1 to another array's row? Let's check.

The line:

aug(i%, 0 TO N%-1) = input(i%,)

This would copy the entire row from input(), which has N% elements (0 to N%-1) into the left half of aug(). Similarly, for the identity part.

But in BBC BASIC, if the input is a 2D array, then input(i%,) refers to the entire row, which is from column 0 to N%-1, which matches the left half. So this line should work.

Then, when creating the identity matrix, setting aug(i%, i% + N%) = 1 is correct.

Then, during normalization:

temp = aug(k%, k%)

aug(k%, k% TO 2*N%-1) = aug(k%, k% TO 2*N%-1) / temp

Yes, because the elements from k% to 2*N%-1 in the pivot row are divided by the pivot element.

For elimination:

For each row i% not equal to k%, compute factor = aug(i%,k%) and subtract factor * pivot row from row i%, starting at column k%.

This should eliminate the element in column k% of row i%.

Using array slices here should be much faster than looping through each element.

Another thing: swapping rows with SWAP. In BBC BASIC, the SWAP statement can swap two variables. But can we swap two array rows with SWAP aug(k%,), aug(pivot_row%,)? I think in BBC BASIC, this is allowed if the arrays are dynamic. Since aug() is a two-dimensional array, but each row is a vector. So SWAP should work for swapping the entire rows.

This would be much faster than looping through each element.

Similarly, when copying the right half of aug() to the output array, the code uses:

output(i%,) = aug(i%, N% TO 2*N%-1)

This copies the slice N% TO 2*N%-1 (which is N% elements) into the output row. Since the output array is N% x N%, this should be correct.

This approach should be efficient as it uses array operations wherever possible, reducing the number of nested loops in BASIC code, which are slower.

Now, testing for possible issues:

What if N% is 0? The code would DIM aug(-1, ...), but the user said "arbitrary-size square matrix", so probably N% >=1. But if N% is zero, the code would error, which is acceptable.

Another thing: the input and output arrays can be the same. In BBC BASIC, if the input and output arrays are the same, then during the copying phase, overwriting the input array with the inverse. But since the code first copies the input into the aug matrix, this is okay. The original input is preserved in the aug matrix until the end, so even if input and output are the same, it's handled correctly.

So this code should work.

Another optimization: instead of using FOR loops to copy the input into the aug matrix and the identity, but using array slices. For example:

aug(, 0 TO N%-1) = input()

But in BBC BASIC, can you assign a block of an array like that? Probably not, unless the left-hand side is a compatible array slice. I think the code as written, using FOR loops to copy each row, is the way to go. Alternatively, if BBC BASIC allows:

FOR i% = 0 TO N%-1

aug(i%, 0 TO N%-1) = input(i%,)

NEXT

Which is what the code does, which is correct.

So overall, this code should correctly invert the matrix using Gaussian elimination with partial pivoting, using array operations where possible for speed.

Here's a BBC BASIC procedure to invert a square matrix using Gaussian elimination with partial pivoting, optimized for speed by utilizing array operations:

bbcbasic
Copy
DEF PROCinvert(input(), output())
LOCAL N%, aug(), i%, k%, pivot_row%, max_val, temp, factor
N% = DIM(input(), 0) + 1

IF DIM(input(), 1) + 1 <> N% THEN ERROR 100, "Matrix not square"

DIM aug(N%-1, 2*N%-1)

FOR i% = 0 TO N%-1
aug(i%, 0 TO N%-1) = input(i%,)
aug(i%, N% TO 2*N%-1) = 0
aug(i%, i%+N%) = 1
NEXT

FOR k% = 0 TO N%-1
pivot_row% = k%
max_val = ABS(aug(k%,k%))
FOR i% = k%+1 TO N%-1
IF ABS(aug(i%,k%)) > max_val THEN
pivot_row% = i% : max_val = ABS(aug(i%,k%))
ENDIF
NEXT

IF max_val = 0 THEN ERROR 101, "Matrix singular"

IF pivot_row% <> k% THEN SWAP aug(k%,), aug(pivot_row%,)

temp = aug(k%,k%)
aug(k%,k% TO 2*N%-1) = aug(k%,k% TO 2*N%-1)/temp

FOR i% = 0 TO N%-1
IF i% <> k% THEN
factor = aug(i%,k%)
aug(i%,k% TO 2*N%-1) -= factor * aug(k%,k% TO 2*N%-1)
ENDIF
NEXT
NEXT

FOR i% = 0 TO N%-1
output(i%,) = aug(i%,N% TO 2*N%-1)
NEXT
ENDPROC
This implementation:

Uses array slicing operations for efficient row manipulation

Performs partial pivoting to ensure numerical stability

Uses in-place operations to minimize memory usage

Processes multiple columns simultaneously using vector operations

Maintains zero-based indexing throughout

The procedure works by:

Creating an augmented matrix with the identity matrix

Performing row operations to convert the left side to identity

Returning the right side of the augmented matrix as the inverse

To use:

Ensure input and output arrays are dimensioned to size [N-1,N-1]

Call with PROCinvert(input(), output())

Handles singular matrices with an error (replace with your own error handling if needed)
Richard Russell
Posts: 272
Joined: Tue 18 Jun 2024, 09:32

Re: Matrix Inversion by DeepSeek

Post by Richard Russell »

Richard Russell wrote: Sat 15 Feb 2025, 16:36 it seemed to think BBC BASIC has array-slicing functionality; I wish that it did but...).
I should have added that BBCSDL and BBCTTY (but not BB4W) do have an array-slicing capability, but implemented in the arraylib.bbc library rather than built-in to the interpreter. Currently that allows only slicing a one-dimensional array but in principle there is no reason why it should not be able to return a whole row, or part of a row, of a two-dimensional array (the slice would still be one-dimensional of course).

I will try to remember to add row-slicing to a future release of the arraylib library. I would suggest a syntax of:

Code: Select all

PROC_rowslice(A(), B(), row%, first%, num%)
where one-dimensional array B() is aliased to part (or all) of a row of two-dimensional array A().