Use Hardware-Efficient Algorithm to Solve Systems of Complex-Valued Linear Equations
This example shows how to solve systems of complex-valued linear equations using hardware-efficient MATLAB® code in Simulink® using a systolic array. The model used in this example shows a hardware-efficient method of solving the system of simultaneous equations
where A is an m-by-n complex matrix, X is an n-by-p complex matrix, and B is an m-by-p complex matrix.
model = 'fi_complex_mldivide_systolic_array_model'; open_system(model)
The enabled subsystems that send the data contain MATLAB Function blocks that keep track of which input to send next, and send the rows of matrices A and B to block 4x4 Complex CORDIC Matrix Left Divide when the
ready signal is high. If you send data when
ready is high, you will not invalidate any data already in the pipeline.
The algorithm overwrites matrix A with an upper-triangular matrix R. The algorithm overwrites B with C = Q'B where Q is unitary and QR = A. The algorithm uses back-substitution on the upper-triangular matrix equation
To examine the algorithm, look Under the mask of block 4x4 Complex CORDIC Matrix Left Divide.
Define complex matrices A and B in the base workspace. In this example, matrix A must be 4-by-4, and matrix B must be 4-by-p, where p is the number of right-hand sides.
rng('default'); A = complex(randn(4,4),randn(4,4)); B = complex(randn(4,1),randn(4,1));
The method uses the CORDIC algorithm, so you must also specify the number of iterations of the CORDIC kernel in the
NumberOfCORDICIterations parameter, or on the block parameters of the 4x4 Complex CORDIC Matrix Left Divide block.
When A and B are double-precision floating-point data types, set the number of CORDIC iterations to the number of bits in the mantissa of a double. If the inputs are fixed point, then the number of CORDIC iterations must be less than the word length. The accuracy of the computation improves one bit for each iteration, up to the word length of the inputs. This model will work with fixed-point, double, and single data types.
NumberOfCORDICIterations = 52;
You must also instantiate variable
BackSubstitutePrototype to specify the data-type used in back-substitution, or enter a prototype value on the block parameters of the "4x4 Complex CORDIC Matrix Left Divide" block. In this case, matrices A and B are double, so set
BackSubstitutePrototype value to be a double. This variable is used by the
cast function to cast the back-substitute variables using the 'like' syntax.
BackSubstitutePrototype = 0;
Turn off expected warnings and simulate the model.
warning_state = warning('off','Coder:builtins:ConstantFoldingOverFlow'); sim(model)
After simulation, the model returns matrix X, which is the solution to the matrix equation
Verify the results by checking that AX-B is a small value.
err = norm(A*X - B)
err = 6.5534e-15
The CORDIC algorithm is used to compute the QR decomposition as in the examples Perform QR Factorization Using CORDIC, and Implement Hardware-Efficient QR Decomposition Using CORDIC in a Systolic Array.
This example is of a 4-by-4 matrix A. You can tile the blocks inside this model to build up to any size matrix.
To see the QR algorithm, look under the mask for block Q'B, R 4x4 Complex CORDIC Systolic Array.
Because the data in this example is complex, an additional step is required to zero out the imaginary parts of the pivots to make them real-valued. This is done by splitting the data into real and imaginary parts and using CORDIC to zero out the imaginary part as if it were two rows of a real matrix. This is equivalent to the complex multiplication
except that the computation is done using the CORDIC algorithm without squares or square-roots.
To see the algorithm for , look under the mask for block Rotate first element to real.
Finally, to compute X, compute the reciprocals of the diagonal elements of R and back-substitute into the right-hand side C. The algorithm uses a CORDIC divide implementation to compute the reciprocals.
To see the back-substitute algorithm, look under the mask for block Back Substitute 4x4 Complex CORDIC.
It is usually unnecessary to explicitly compute the inverse of a matrix ,. However, if you want to do so, set B equal to the identity matrix I. Then, the solution of the equation
A = complex(2*rand(4,4)-1,2*rand(4,4)-1); B = complex(eye(4)); NumberOfCORDICIterations = 52; BackSubstitutePrototype = 0;
Simulate the model
Verify that AX and XA are close to the identity matrix. There will be small differences due to round-off errors.
ans = 1.0000 + 0.0000i 0.0000 - 0.0000i 0.0000 + 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i 1.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i 0.0000 - 0.0000i 1.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i 1.0000 - 0.0000i
ans = 1.0000 - 0.0000i -0.0000 - 0.0000i -0.0000 + 0.0000i 0.0000 - 0.0000i 0.0000 + 0.0000i 1.0000 + 0.0000i 0.0000 - 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 1.0000 - 0.0000i -0.0000 + 0.0000i 0.0000 - 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i 1.0000 + 0.0000i
Normalizing to fractional types makes the computations easier in fixed-point. You can scale the matrices so that their data is in the range [-1, +1] and use fractional fixed-point types because the solution to the matrix equation
is the same as the solution to
To convert to fractional scaling with no additional cost in generated code, you can use the
reinterpretcast function in MATLAB or the Data Type Conversion block in Simulink with Input and output to have equal Stored Integer.
Run Exhaustive Test Points
You can run many test inputs through the model by making A and B three-dimensional arrays.
m = 4; n = 4; p = 1; n_test_inputs = 100;
Create the inputs such that the real and imaginary parts are in the range [-1, +1]
A = complex(2*rand(m,n,n_test_inputs)-1, 2*rand(m,n,n_test_inputs)-1); B = complex(2*rand(m,p,n_test_inputs)-1, 2*rand(m,p,n_test_inputs)-1);
In this example, matrices A and B contained real and imaginary parts already in the range [-1, +1], so set the types to be fractional.
data_word_length = 24; data_fraction_length = 23;
QR Data Types
The growth in the elements of R in the real QR factorization is (see Perform QR Factorization Using CORDIC). The elements are complex in this example, so an additional growth factor of is needed because
Also, the CORDIC algorithm grows intermediate values by the following gain factor where is the number of CORDIC iterations before it is normalized out.
Therefore, an upperbound for growth in the QR algorithm for m-by-n complex matrix A is the product of all the growth factors:
It this example, A is 4-by-4, so the maximum growth factor in the QR algorithm is
Therefore, the number of additional integer bits to allow in the QR algorithm to avoid overflow when you have m=4 rows is:
qr_growth_bits = ceil(log2(1.6468*sqrt(2*m)))
qr_growth_bits = 3
Cast Matrices A and B to Fixed Point
The model requires that the inputs be signed and the word length of B be the same as the word length of A.
Grow the data wordlength to accommodate the QR growth.
qr_input_word_length = data_word_length + qr_growth_bits
qr_input_word_length = 27
Cast A to fixed point and B to A's type.
T.A = fi(, 1, qr_input_word_length, data_fraction_length); T.B = T.A; A = cast(A,'like',T.A); B = cast(B,'like',T.B);
Back-Substitute Data Type
If is an -by- invertible complex matrix, and is the solution of matrix equation , then using properties of vector and matrix norms (see, for example, reference ), it can be shown that
where and is the smallest singular value of . This bound is related to the condition number of .
In this example, is a complex vector with the real and imaginary parts bounded by 1. Therefore, .
Therefore, if you know the distribution of singular values for the matrices in your problem, then you can choose the number of additional integer bits required to avoid overflow in the back-substitute with a given probability (see, for example, reference ).
In this example, we compute the singular values for all the matrices in our test bench and choose the number of integer bits to avoid overflow based on that.
singular_values = zeros(n,n_test_inputs); A0 = double(A); for k = 1:n_test_inputs singular_values(:,k) = svd(A0(:,:,k)); end condition_numbers = singular_values(1,:)./singular_values(end,:); x_bound = sqrt(2*n)./singular_values(end,:);
The number of bits of growth to add is the base-2 logarithm of the maximum bound.
integer_bits_for_backsubstitute = ceil(log2(max(x_bound)))
integer_bits_for_backsubstitute = 7
Subtract the required integer bits for back-substitute from the wordlength, and an additional bit for the sign bit.
backsubstitute_fraction_length = T.A.WordLength - integer_bits_for_backsubstitute - 1
backsubstitute_fraction_length = 19
Therefore, the data type for the back-substitute is:
BackSubstitutePrototype = fi(0, 1, T.A.WordLength, backsubstitute_fraction_length)
BackSubstitutePrototype = 0 DataTypeMode: Fixed-point: binary point scaling Signedness: Signed WordLength: 27 FractionLength: 19
Set the number of CORDIC iterations to be one less than the fixed-point word length of A.
NumberOfCORDICIterations = T.A.WordLength - 1
NumberOfCORDICIterations = 26
Simulate the model with fixed-point inputs.
Calculate and Plot the Errors
A measure of error is
norm(AX - B)
for each pair of inputs A and B.
norm_error = zeros(1,size(X,3)); B0 = double(B); X0 = double(X); for k = 1:size(X,3) norm_error(k) = norm(A0(:,:,k)*X0(:,:,k) - B0(:,:,k)); end
Plot the errors. The errors are low because of the data types chosen. The errors are typically higher when the condition number of matrix A is high, as the theory predicts.
figure(1) clf h1 = subplot(2,1,1); plot(norm_error) grid on title('Errors') ylabel('norm(AX - B)') h2 = subplot(2,1,2); plot(condition_numbers) grid on title('Condition numbers') ylabel('cond(A)') xlabel('Test points') linkaxes([h1,h2],'x')
 Zizhong Chen and Jack J. Dongarra. "Condition Numbers of Gaussian Random Matrices". SIAM Journal on Matrix Analysis and Applications. 27.3 (July 2005), pp. 603-620.
 George E. Forsythe and Cleve B. Moler. Computer Solution of Linear Algebraic Systems. Englewood Cliffs, N.J.: Prentice-Hall, 1967.
 George E. Forsythe, M.A. Malcom and Cleve B. Moler. Computer Methods for Mathematical Computations. Englewood Cliffs, N.J.: Prentice-Hall, 1977.
 Cleve B. Moler. Cleve's Corner: What is the Condition Number of a Matrix?>, The MathWorks, Inc. 2017.
 Cleve B. Moler. Numerical Computing with MATLAB. SIAM, 2004. isbn: 978-0-898716-60-3.
 Gene H. Golub and Charles F. Van Loan. Matrix Computations. The Johns Hopkins University Press.