key: cord-0059028-uqwpnaqw
authors: Skala, Vaclav
title: Conditionality of Linear Systems of Equation and Matrices Using Projective Geometric Algebra
date: 2020-08-19
journal: Computational Science and Its Applications - ICCSA 2020
DOI: 10.1007/978-3-030-58802-1_1
sha: 9957aec4845861fb38882242b7af4f152fe71a89
doc_id: 59028
cord_uid: uqwpnaqw

Linear systems of equations and their reliable solution is a key part of nearly all computational systems and in a solution of many engineering problems. Mostly, the estimation of the matrix conditionality is used for an assessment of the solvability of linear systems, which are important for interpolation, approximation, and solution of partial differential equations especially for large data sets with large range of values. In this contribution, a new approach to the matrix conditionality and the solvability of the linear systems of equations is presented. It is based on the application of the geometric algebra with the projective space representation using homogeneous coordinates representation. However, the physical floating-point representation, mostly the IEEE 754-219, has to be strictly respected as many algorithms fail due to this.

Solutions of linear systems of equations are probably one of the most used operations in many applications. Generally, there are two main groups of those:

-non-homogeneous systems of linear equations, i.e. Ax = b -homogeneous system of equations, i.e. Ax = 0

Those two groups of linear systems are solved differently, in spite of the fact, that they are dual problems in some sense. Therefore, there should be common methods for solving both groups of such linear systems of equations. The first case is deeply analyzed in the courses of numerical mathematics, while the second one is considered as marginal without advanced solution methods. (In the following, standard properties of the matrices are expected).

Using the principle of duality and projective extension of the Euclidean space the first type of the linear system, i.e. Ax = b, can be easily transformed to the second type, i.e. Ax = 0. The geometric algebra offers more general formalism, which can be used for a better understanding of the linear system of equations properties and behavior.

However, in practical applications, any implementation has to take into consideration the IEEE 754-219 or similar standards of the floating-point physical representation.

Both types of the linear systems of equations, i.e. Ax = b and Ax = 0, have a simple geometrical interpretation, which can be easily demonstrated on two simple problems in the E 2 case, see Fig. 1 :

-an intersection computation of two lines, i.e. Ax = b; the matrix A is n × n -a line determination given by two points, i.e. Ax = 0; the matrix A is n × (n + 1)

Similarly, in the E 3 case, where three planes intersect in a point or three points define a plane, if no singular cases are considered. Also, the singular or close to the singular case has to be properly solved in both cases. In the first case, if the lines are collinear, in the second one, if two points are the same.

Unfortunately, programmers do not solve such cases correctly in many cases. Usually, a condition like det(A) <= eps is used for a "singularity" detection and eps is taken as a "reasonable" small number.

Let us consider the intersection computation of two lines in the E 2 for simplicity. The Cramer's rule can be used and the solution of Ax = b is given as x = Det x /Det and y = Det y /Det. If Det x ∼ det(A) then the x value might be close to 1 regardless of the eps value.

As those two problems in Fig. 1 are dual, i.e. the meet operator in the first case and union operator in the second one, the principle of duality must be kept [7] . Similarly, in the E 3 case in Fig. 2 . 

The vector algebra (Gibbs algebra), which is used in engineering practice, uses two basic operations on two vectors a, b in E n :

-the inner product (scalar product or dot product) c = a · b, where c is a scalar value -the outer product c = a ∧ b (the cross-product in E 3 c = a × b), where c is a bivector and has different properties as it represents an oriented area in the n-dimensional space, in general.

The Geometric Algebra (GA) uses a "new" product called geometric product defined as:

where ab is the new entity. It should be noted, that it is a "bundle" of objects with different dimensionalities and properties, in general.

In the case of the n-dimensional space, the vectors are defined as a = (a 1 e 1 + ... + a n e n ), b = (b 1 e 1 + ... + b n e n ) and the e i vectors form orthonormal vector basis in E n . In the E 3 case, the following objects can be used if geometric algebra [26] :

The significant advantage of the geometric algebra is, that it is more general that than the Gibbs algebra and can handle all objects with dimensionality up to n.

The geometry algebra uses the following operations, including the inverse of a vector.

It should be noted, that geometric algebra is anti-commutative and the "pseudoscalar" I in the E 3 case has the basis e 1 e 2 e 3 (briefly as e 123 ), i.e.

e i e j = −e j e i e i e i = 1 e 1 e 2 e 3 = I a ∧ b ∧ c = q

where q is a scalar value (actually a pseudoscalar).

The geometric product of two vectors is represented as:

It is not a "friendly user" notation for a practical application and causes problems in practical implementations, especially due to anti-commutativity of the geometric product. However, the geometric product can be easily represented by the tensor product, which can be easily represented by a matrix. The tensor product for the 4-dimensional case, as the homogeneous coordinates will be used in the following, is represented as:

where B + U + D are Bottom triangular, Upper triangular, Diagonal matrices, a 4 , b 4 will be used for the homogeneous coordinates, i.e. actually w a , w b (the projective notation will be explained later), and the operator ⊗ means the anti-commutative tensor product.

Let us consider the projective extension of the Euclidean space in the E 3 case and use of the homogeneous coordinates. The "hidden" advantage of this is that it enables us to represents points or values close to infinity. Let us consider a vector a = [a 1 , a 2 , a 3 : a 4 ] T , a 4 = 0, which represents actually the vector (a 1 /a 4 , a 2 /a 4 , a 3 /a 4 ) in the E 3 space. The ":" symbol within the projective representation is used to distinguish different properties. In the case of a plane, the vector (a 1 , a 2 , a 3 ) represents the "normal vector" of a plane ( actually it is a bivector), while a 4 actually represents a pseudo-distance of a plane from the origin.

Let us consider a second vector

It can be seen, that the diagonal of the matrix Q actually represents the inner product in the projective representation as:

where " ≡ " means projectively equivalent. In the E 3 case, the inner product actually projectively represents the trace tr(Q) of the matrix Q. The outer product (the crossproduct in the E 3 case) is then represented respecting the anti-commutativity as:

It should be noted, that the outer product can be used for a solution of a linear system of equations Ax = b or Ax = 0, too.

The principle of duality is an important principle, in general. The duality principle for basic geometric entities and operators is presented by Table 1 and It means, that in the E 2 case a point is dual to a line and vice versa, the intersection of two lines is dual to a union of two points, i.e. line given by two points; similarly for the E 3 case, where a point is dual to a plane [7, 14, 15] . The direct consequence of the principle of duality is that, the intersection point x of two lines p 1 , p 2 , resp. a line p passing two given points x 1 , x 2 , is given as:

similarly in the dual case.

In the case of the E 3 space, a point is dual to a plane and vice versa. It means that the intersection point x of three planes ρ 1 , ρ 2 , ρ 3 , resp. a plane ρ passing three given points x 1 , x 2 , x 3 is given as:

It can be seen that the above formulae is equivalent to the "extended" cross-product, which is natively supported by the GPU architecture. For an intersection computation, we get:

Due to the principle of duality, the dual problem solution is given as:

The above presented formulae proved the strength of the formal notation of the geometric algebra approach. Therefore, there is a natural question, what is the more convenient computation of the geometric product, as computation with the outer product, i.e. extended cross product, using standard vector algebra approach is not simple.

Fortunately, the geometric product of ρ 1 , ρ 2 , resp. of x 1 and x 2 vectors using homogeneous coordinates given as the anti-commutative tensor product is given as:

The question is how to compute a line p ∈ E 3 given as an intersection of two planes ρ 1 , ρ 2 , which is the dual problem to a line given by two points x 1 , x 2 as those problems are dual. It leads to the Plücker coordinates used especially in robotics [17, 21] . Using the geometric algebra, the principle of duality and projective representation, we can directly write:

It can be seen that the formula given above keeps the duality in the final formulae, too. A more detailed description can be found in [17] . It should be noted that the geometric algebra and the principle of duality offers:

-to solve dual problems by one programming sequence -it naturally supports parallel programming as it uses vector-vector operations as the SSE instructions or GPU can be used -the solution can avoid a division operation if the result can be left in the projective representation [24] -results of operations are in symbolic form, which enables further symbolic processing using vector algebra or geometric algebra rules

A solution of a linear system of equations is a part of the linear algebra and used in many computational systems. There are many publications related to methods for a linear system of equations following different criteria:

-memory requirements -linear systems with sparse or tri-diagonal matrices [8] , etc.

-algorithm computational complexity [25] using a block decomposition -parallel processing of large linear systems, especially in connection with GPUs applications [28] -higher precision of computation and numerical stability [3, 24] .

It should be noted, that (ξ 1 , ..., ξ n ) gives the "direction", only.

As the solution of a linear system of equations is equivalent to the outer product (generalized cross-vector) of vectors formed by rows of the matrix D, the solution of the system Dξ = 0 is defined as:

where:

The application of the projective extension of the Euclidean space enables us to transform the non-homogeneous system of linear equations Ax = b to the homogeneous linear system Dξ = 0, i.e.:

It is a very important result as a solution of a linear system of equations is formally the same for both types, i.e. homogeneous linear systems Ax = 0 and non-homogeneous systems Ax = b. It means, that many problems might be transformed to the outer product application, e.g.. computation of barycentric coordinates [13] , length, area, volume og geometric entities [12] , etc.

As the solution is formally determined by Eq. 14, formal linear operators can be used for further symbolic processing using standard formula manipulation, as the geometry algebra is multi-linear. Even more, it is capable to handle more complex objects generally in the n-dimensional space, i.e. oriented surfaces, volumes, etc.

Also, it is possible to use the 'functional analysis approach: "Let L is a linear operator, then the following operation is valid....". As there are many linear operators like derivation, integration, Laplace transform, etc., there is a huge potential of applications of that to the formal solution of the linear system of equations.

However, it is necessary to respect, that in the case of projective representation specific care is to be taken for deriving the rules for derivation, etc., as actually, a fraction is to be processed.

It should be noted, that a solution is in symbolic representation, too, which can be used for further processing and symbolic manipulation.

Both types of the linear systems of equations, i.e. Ax = b (A is n × n) and Ax = 0 (A is (n+1)×n), actually have the same form Ax = 0 (A is (n+1)×n), now, if the projective representation is used. Therefore, it is possible to show the differences between the matrix conditionality and conditionality (solvability) of a linear system of equations, see Fig. 3 .

For the matrix conditionality, the eigenvalues are usually used [10, 11] and the ratio rat λ = |λ max |/|λ min | is mostly used as a criterion. It should be noted, that λ ∈ C, i.e the eigenvalues are complex, in general. If the ration rat λ is high, the matrix is said to be ill-conditioned [4, 6] , especially in the case of large data with a large span of data [22, 23] . There are two cases, which are needed to be taken into consideration:

-non-homogeneous systems of linear equations, i.e. Ax = b. In this case, the matrix conditionality is considered as a criterion for the solvability of the linear system of equations. It depends on the matrix A properties, i.e. on eigenvalues. 

A conditionality number κ(A) = |λ max |/|λ min | is usually used as the solvability criterion. In the case of the Eq. 16, the matrix conditionality is κ(A) = 10 2 /10 −2 = 10 4 . However, if the 1 st row is multiplied by 10 −2 and the 3 rd row is multiplied by 10 2 , then the conditionality is κ(A) = 1. Therefore, the estimation of the matrix conditionality using ratio of eigenvalues is not applicable. It should be noted, that the right side vector b is to be appropriately multiplied. Classification of possible conditionality evaluation:

• volumetric based -as a matrix conditionality the value of det(A) is taken. However, its value is given by a volume of a hyper-polytope (n-dimensional parallelepiped) given by the rows of the matrix A. The "standard" matrix conditionality using eigenvectors is not quite appropriate, as shown above, Eq. 16. • area based -the area of bivectors defined as ||a i ∧ a j ||. This evaluation was recently used for higher size of the Hilbert matrix is [19, 20] . However, the results were not satisfactory in some cases. 

For a solution of a linear system of equations some specific operations are used, e.g.. row multiplication by a non-zero constant, which does not have any influence to the result correctness. Let us have a look at such operations more deeply. The basic operations can be split into two types:

row multiplication by a constant p = 0. It changes the row values, but does not change the result of the linear system Ax = 0, resp. Ax = b, if the vector b is multiplied appropriately. column multiplication -by a constant s = 0. This operation is not usually mentioned in numerical methods. It means, that there is a scaling of the relevant element of the vector x of the system Ax = 0, resp. Ax = b.

It means, that those two main operations with both linear system of equations Ax = 0, resp. Ax = b can be described as

resp.

where: P, resp. S are diagonal matrices with p i = 0, s j = 0 multiplicative constants. The matrix S represents an-isotropic scaling of the size (n + 1) × (n + 1) in the first case, resp. n × n in the second one. It should be noted, that the diagonal matrix P does not have any influence on the result, due to the projective representation, and because S S −1 is the identity matrix the solution of the linear system of equations is correct in both cases. It can be easily proved that both operations do have an influence on eigenvalues and therefore κ(A) = |λ max |/|λ min | cannot be used for assessment of the solvability of the linear system of equations. For the same reasons, also the assessment based of the area on bivectors cannot be used.

The only angular criterion is invariant to the row multiplications, while only the column multiplication changes angles of the bivectors. There are several significant consequences for the numerical solution of linear systems of equations:

-the solvability of a linear system of equations can be improved by the column multiplications, only, if unlimited precision is considered. Therefore, the matrix-based pre-conditioners might not solve the solvability problems and might introduce additional numerical problems. -the precision of computation is significantly influenced by addition and subtraction operations, as the exponents must be the same for those operations with mantissa. Also, the multiplication and division operations using exponent change by 2 ±k should be preferred.

However, the question is, whether the above presented operations cannot be efficiently used for better numerical stability of numerical solutions. In the matrix operations, it is necessary to respect the floating-point representation [27] as it is a crucial factor in the solution of large or ill-conditioned linear systems of equations.

There are several methods used to improve the ratio κ(A) = |λ max |/|λ min | of the matrix A of the linear system, e.g.. matrix eigenvalues shifting or preconditioning [1, 2] . The preconditioning is usually based on solving a linear system Ax = 0:

where M −1 is a matrix, which can cover complicated computation, including Fourier transform. The inverse operation, i.e. M −1 , is computationally very expensive as it is of O(n 3 ) complexity. Therefore, they are not easily applicable for large systems of linear equations used nowadays. There are some methods based on incomplete factorization, etc., which might be used [11] .

The proposed matrix conditionality improvement method requires only determination of the diagonal matrices values P and S, i.e. multiplicative coefficients p i = 0, s j = 0, which have to be optimized. This is a significant reduction of computational complexity, as it decreases the cost of finding sub-optimal p i , s j values.

The proposed approach was tested on the Hilbert's matrix, which is extremely illconditioned. The Hilbert matrix conditionality can be estimated as [5, 9] .

κ(H n ) e 3.5n (20) where n is the size of the Hilbert matrix; for n = 10 the conditionality is estimated as κ(H n ) 1.56810 15 .

The Hilbert matrix H defined as:

where: i, j = 1, ...n.

The experimental results of the original conditionality κ(H orig ) and conditionality using the proposed method κ(H new ) are presented in Table 3 .

To prove the proposed approach, the Hilbert matrix was taken, as it is ill-conditioned. The conditionality of the matrix H orig , modified by diagonal matrices P and S, was evaluated by the cond(matrix) function in the Octave system using the simplified Monte Carlo approach. The experiments proved, that the conditionality cond(H new ) of the modified matrix using the proposed approach was decreased by more than half of the magnitude for higher values of n, see Table 3 . This is consistent with the recently obtained results [16] , where the inverse Hilbert matrix computation using the modified Gauss elimination without division operation was analyzed.

The Hilbert matrix conditionality improvement also improved the angular criterion based on maximizing the ratio κ rat (H) defined as: κ rat (H) = cos γ min cosγ max (22) It says, how the angles cos γ i j , formed by the vectors a i j of the bivectors are similar, see Fig. 3 . It means, that if the ratio κ rat (A) 1 the angles of all bivectors are nearly equal.

In the case of conditionality assessment of the linear system of equations Ax = 0, the angles β i j , formed by the angels α i j have to be taken into account, see Fig. 3 . The ratio κ rat (H) is then defined as:

The advantage of the angular criterion is that it is common for the conditionality evaluation in the cases, i.e.:

-the matrix conditionality -conditionality of the linear system of equations It should be noted, that this conditionality assessment method gives different values of conditionality of those two different cases, as in the first case only the matrix is evaluated, while in the second one the value of the b in the Ax = b is taken into account. The results presented in Table 4 reflects the improvement of the Hilbert matrix by proposed approach using the diagonal matrices P and S used as the multipliers.

This contribution briefly presents a new approach to the matrix conditionality and the solvability of linear systems of equations assessment, based on the bivectors angular properties based on geometry algebra. The presented approach enables to make conditionality assessment for both types of linear equations, i.e. Ax = b and Ax = 0.

Also, the equivalence of the solution of linear systems of equations with the application of the outer product (extended cross product) has been presented. It offers simple and efficient solutions to many computational problems, if combined with the principle of duality and projective notation. Also, the presented approach supports direct GPU application with a potential of significant speed-up and parallelism. Even more, the approach is applicable to n-dimensional problem solutions, as the geometric algebra is multidimensional.

In the future, the proposed approach will be more deeply analyzed and experimentally verified using numerical operations with values restricted to 2 k as in this case the mantissa value is not changed. The expected application is in the radial basis function interpolation and approximation for large data, when matrices are ill-conditioned, and in the solution of partial differential equations. Now, reduction can be used ||a i ∧ a j || 2 = ||a i || 2 ||a j || 2 − (a i · a j ) 2 (8) and finally the square of the area of a bivector is given as:

Preconditioning techniques for large linear systems: a survey

Matrix Preconditioning Techniques and Applications. Cambridge Monographs on Applied and Computational Mathematics

Preconditioning for sparse linear systems at the dawn of the 21st century: history, current developments, and future perspectives

The conditionally and expected error of two methods of computing the pseudo-eigenvalues of a complex matrix

Accuracy and Stability of Numerical Algorithms

On the conditionality of eigenvalues close to the boundary of the numerical range of a matrix

Proof by duality: or the discovery of "new

Big geodata surface approximation using radial basis functions: a comparative study

Scientific Computing with MATLAB and Octave

Iterative Methods for Sparse Linear Systems

Numerical problems for Large Eigenvalue Problems

Length, area and volume computation in homogeneous coordinates

Barycentric coordinates computation in homogeneous coordinates

Intersection computation in projective space using homogeneous coordinates

Duality, barycentric coordinates and intersection computation in projective space with GPU support

Modified Gaussian elimination without division operations

Plücker Coordinates and Extended Cross Product for Robust and Fast Intersection Computation

Extended cross-product" and solution of a linear system of equations

High dimensional and large span data least square error: numerical stability and conditionality

Least square method robustness of computations: what is not usually considered and taught

Projective Geometry, Duality and Plücker Coordinates for Geometric Computations with Determinants on GPUs

RBF Interpolation with CSRBF of Large Data Sets

RBF Approximation of Big Data Sets with Large Span of Data

A Precision of Computation in the Projective Space

Gaussian elimination is not optimal

Geometric Algebra for Computer Graphics

Distributed-memory lattice h-matrix factorization

Acknowledgment. The author would like to thank their colleagues and students at the University of West Bohemia, Plzen, for their discussions and suggestions and implementations, especially to Michal Smolik and Martin Cervenka for discussions and help with the Octave implementation, to Filip Hacha, Vitek Poor and Jan Kasak for additional verification of some parts in the Octave system. Thanks belong also to anonymous reviewers for their valuable comments and hints provided. This research was supported by the Czech Science Foundation (GACR) No. GA 17-05534S.

Computation of a bivector area in the n-dimensional space can be made as follows:Therefore the square of the bivector area is given asAs the following identity is validthen it can be expressed asby a substitution of Eq. 1and algebraic manipulationi.e.||a i ∧ a j || 2 = ||a i || 2 ||a j || 2 ||a i || 2 ||a j || 2 − (a i · a j ) 2 ||a i || 2 ||a j || 2