A Generalized Matrix Inverse with Applications to Robotic Systems

05/07/2018 ∙ by Bo Zhang, et al. ∙ 0

It is well-understood that the robustness of mechanical and robotic control systems depends critically on minimizing sensitivity to arbitrary application-specific details whenever possible. For example, if a system is defined and performs well in one particular Euclidean coordinate frame then it should be expected to perform identically if that coordinate frame is arbitrarily rotated or scaled. Similarly, the performance of the system should not be affected if its key parameters are all consistently defined in metric units or in imperial units. In this paper we show that a recently introduced generalized matrix inverse permits performance consistency to be rigorously guaranteed in control systems that require solutions to underdetermined and/or overdetermined systems of equations.



There are no comments yet.


page 8

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

I Introduction

Many robotic and mechatronic control systems are mathematically represented, analyzed, and ultimately implemented as compositions of linear systems. This is the case even if the dynamics of the system are fundamentally nonlinear but are solved in terms of linear-algebraic equations or locally-linear approximations within globally nonlinear state-space models. In such systems it is commonly necessary to solve an overdetermined or underdetermined set of equations in order to satisfy a given set of constraints or to select from a multiplicity of local solutions within an iterative process. Although this may happen in a single step of a logically small component of a very large system, if care is not taken to ensure that critical mathematical properties are properly preserved the integrity of the overall system can be compromised.

It is a fundamental design principal that sensitivity to arbitrary application-specific details should be minimized whenever possible. For example, if a system is defined and performs well in some particular Euclidean coordinate frame then it should be expected to perform identically if that coordinate frame is arbitrarily rotated or scaled. Similarly, the performance of the system should not be affected if its key parameters are all consistently defined in metric units or in imperial units. In this paper we show that a recently introduced generalized matrix inverse permits performance consistency to be rigorously guaranteed in control systems that require solutions to underdetermined and/or overdetermined systems of equations. Specifically, we show that consistency with respect to arbitrary choices of units in state-space models of robotic systems can be affected by a simple replacement of the Moore-Penrose generalized matrix inverse with a general unit-consistent inverse.

Ii Generalized Matrix Inverses

For a nonsingular matrix there exists a unique matrix inverse, , which preserves many properties that hold for ordinary scalar inverses, e.g., matrix inversion distributes over nonsingular multiplicands as:


where noncommutativity of matrix multiplication imposes a constraint on the ordering of terms but is otherwise analogous to the scalar case. In a practical application the above inverse-distributivity property implies that if we only have access to

from its inverse in a linearly transformed space,

, then can be obtained simply as .

When attempting to generalize the notion of a matrix inverse for singular it is only possible to define an approximate inverse that retains a subset of the algebraic properties of a true matrix inverse [1], such as:




and/or other properties that may be of analytic or application-specific utility. The Moore-Penrose pseudoinverse [2, 3] (MP inverse), , is by far the most widely known and used generalized inverse111The Matlab/Octave operator pinv(M) returns the MP inverse of its matrix argument.. It is defined for any matrix and satisfies the above generalized inverse properties as well as the following for any conformant unitary/orthonormal222For notational purposes we retain the generality of interpreting and as arbitrary unitary matrices over or (e.g., is equal to its conjugate transpose ), but for all practical purposes in this paper they can be thought of as representing permutations of state variables and/or rotations of a global coordinate frame. matrices and :


This property implies that the MP inverse is applicable to problems defined in a Euclidean state space for which the behavior of the system of interest should be invariant with respect to arbitrary rotations of the coordinate frame. In that context can be recovered from its MP inverse in a rotationally transformed space as , where consistency with respect to rigid rotations has been implicitly exploited. This can be understood by noting that if is singular it must be assumed that for arbitrary nonsingular matrices and . More technically, the MP inverse is consistent with respect to arbitrary unitary transformations but not to general linear transformations.

Despite its widespread default use throughout most areas of engineering (often implicitly under the name “least-squares”), the MP inverse does not satisfy conditions appropriate for many problems to which it is commonly applied, e.g., ones that require consistency with respect to the choice of units for state variables. For example, a state parameterized with four variables defined respectively in units relating to temperature, pressure, speed, and distance can be thought of as defining a 4-dimensional Cartesian coordinate frame, but it would make no sense to rotate that coordinate frame to a space in which these variables are mixed. In this case consistency should be preserved with respect to changes of units, e.g., from imperial to metric, rather than with respect to rotations of a global coordinate frame which has no physical meaning or interpretation. This kind of unit consistency (UC) requires a generalized inverse that satisfies


where the diagonal matrix represents units on variables in one space and the diagonal matrix represents different units for the same variables in a different space.

The hazards associated with the misuse of the MP inverse have been noted in the robotics literature [4, 5], and disciplined methodologies have been developed to address the issue in common situations that arise in that context [6, 7]. However, the recent derivation of an inherently unit-consistent generalized inverse, or UC inverse, reduces the need for tailored solutions because in principle it can be simply substituted in place of the MP inverse [8]. This is not only useful because it simplifies the implementation process, it also reduces the opportunity for subtle implementation errors to be introduced.

In the next section we examine a robotic system in which we demonstrate that the MP inverse fails to preserve unit consistency and thus produces unreliable results. We then show that simply replacing the MP inverse with the UC inverse provides improved and completely stable behavior that is invariant with respect to arbitrary changes of units on key parameters.

Iii Generalized Matrix Inverse for Robotics System

In this section we consider an example involving the motion of a robotic arm. This example is based on systems that have been studied in the literature to show the important aspects of real-world robotics and mechanical system applications [7, 9]. We will start by describing the structure of the robotic arm and the equations that model its motion so that it can be controlled to perform desired operations.

The components and configuration of the robotic arm are shown in Figure 1. It has structures common in mechanical systems. The tip-point

is the end-effector, which denotes the end of a robotic arm and designed to interact with the outside environment. Each joint can be actuated by a motor, and due to the design of the connection, it is limited to 2-dimensional planar motion. The mechanism has 3 degrees of freedom: two of them are rotational and one is linear. These parameters are initialized as

, and . A mixed control of these joints will determine the motion of tip-point , e.g., to take it to a desired state B. Taking the cylindrical joint on the frame as the origin, the position of the tip-point () is given by equations 6-8.


where . A Jacobian matrix represents the transformation of the end-effector representation of the dynamic system (

) into a joint-state representation in which all of the key parameters are represented as a vector as

. For the current problem the state of the system is given by , which are the states of the 3 joints. The Jacobian matrix for this system,


is determined to be

Fig. 1: Robotic arm system with two rotational joints and one linear joint. The tip-point is the end-effector designed to interact with the environment.

The target velocity is thus




where is the derivative of , and is the joint velocity of the robotic arm represented by . Assuming a target velocity , the solution to achieve the desired motion is


but this cannot be evaluated if is singular, and it will be singular in this case because the motion is constrained to a 2-dimensional plane in a 3-dimensional space. The solution therefore requires use of a generalized matrix inverse in place of the undefined matrix inverse. As has been discussed, the most commonly used generalized inverse is the Moore-Penrose (MP) inverse, which gives a result for as


where is the MP inverse. As will be seen, the use of the MP inverse is not appropriate if the state vector has variables with units that must be preserved, which is true for the case of . In this system lengths are defined in meters, , , , and , though the system should be expected to work correctly no matter what units are used as long as everything is defined consistently in those units.

The joint velocity needed to keep the speed of tip-point equal to is calculated using equation 14 with timesteps of over a simulation time of 0.1 s. At each timestep (e.g., ) the angles and lengths are updated for the next iteration as so that the Jacobian matrix can also be updated to solve the for the next timestep. The calculated velocities included in Table I can be seen to change slowly over time, and the simulation shows that they lead to a tip-point velocity that is approximately equal to during the process. The final state of the tip-point is shown in Figure 2(a), which shows that the target was successfully controlled to move along the designed track.

0.000 -27.881 -12.12 -1.543 [2, -2, 0]
0.001 -27.826 -11.981 -1.548 [2, -2, 0]
0.002 -27.772 -11.838 -1.553 [2, -2, 0]
0.003 -27.719 -11.695 -1.558 [2, -2, 0]
0.004 -27.666 -11.552 -1.563 [2, -2, 0]
0.005 -27.614 -11.409 -1.568 [2, -2, 0]
0.006 -27.563 -11.266 -1.573 [2, -2, 0]
0.007 -27.513 -11.123 -1.578 [2, -2, 0]
0.008 -27.464 -10.980 -1.582 [2, -2, 0]
0.009 -27.414 -10.837 -1.587 [2, -2, 0]
0.010 -27.367 -10.693 -1.592 [2, -2, 0]
TABLE I: Joint velocity of robotic arm calculated using MP inverse. Timestep: 0.001s.

In contrast to the first simulation, the second simulation is performed identically except that lengths are defined in units of centimeters instead of meters: , , , and . As a result, some elements of the Jacobian matrix are changed in magnitude, but the overall behavior of the system should be expected to remain unchanged. In other words, as long as everything is implemented consistently using lengths defined in meters, or implemented consistently using lengths defined in centimeters, the behavior of the system should be the same because the controls will be computed accordingly. The choice of units should not matter, and if it does then the system cannot be trusted.

Although the joint velocities are calculated the same way as before with a timestep of , the system now diverges rapidly. For example, Figure 3(a) shows that the value of the angular velocity rapidly fluctuates between , whereas the computed angular velocity in the previous simulation remained stable at around . Figure 2(b) shows that the control is unsuccessful as the tip-point motion diverges from the designed path. To reduce the magnitude of this deviation the control system has to be performed using much smaller timesteps, which requires control operations to be calculated and applied much more frequently to keep the tip from diverging too far from the correct path.

Fig. 2: Tip point position after simulation, solved with MP inverse. (a) length unit(m), timestep (). (b) length unit(cm), timestep (). (c) length unit(cm), timestep (). Simulation (b) fails by deviating from designed track.

The results of several tests are plotted in Figure 3. It can be seen that acceptable results are obtained with a timestep of when lengths are defined in units of meters, but unstable results are produced when the length unit is changed to centimeters. In the centimeter case it is necessary to reduce the timestep by an order of magnitude to to achieve reliable control, but this increases the computational complexity by an order of magnitude because control operations must be calculated much more frequently. While this does yield an acceptable tip-point velocity, the control values are different from those produced when units were in meters and this is evidence that something is wrong. It is tempting to conclude based on the final state displayed in Figure 2(c) that the controls perform equally well in both cases because they produce seemingly similar paths. However, Figure 3(a) shows that the control produced in the case of centimeter units leads to very erratic high-frequency motion along the path, and the controls required to achieve these wildly changing velocities may be physically difficult to actually achieve with real hardware. Even if they can be achieved they may cause excessive stress and wear on machine elements.

The simulations clearly show that different units lead to different configurations for the same problem when using the MP inverse. The final state of the robotic arm is , and for lengths in meters, , and for lengths in centimeters. The fact that results depend on the choice of units raises concerns because there is no way to predict when and how this dependency may lead to bad results. Although the timestep can be reduced until the controls calculated using the MP inverse produce acceptable results, the need to adjust the timestep to compensate for arbitrary dependencies on the choice of units implies that the solution has been tailored to the problem and the results may not remain acceptable for a slightly different configuration. Therefore, this control system using the MP inverse cannot be trusted.

Fig. 3: Joint velocity () solved with MP inverse. (a) length unit(cm), timestep (). (b) length unit(cm), timestep (). (c) length unit(cm), timestep (). (d) length unit(m), timestep (). Simulation (a) does not converge, and simulation (b) and (c) can converge with a smaller timestep.

It is important to review why the control system in this case seemed to work well when units were defined in meters but not when they were defined in centimeters. This happened because the control calculations used the MP inverse when there was need to invert a singular Jacobian matrix. Based on the properties discussed in the previous chapter, the MP inverse is defined to provide a solution to a linear system that minimizes the Euclidean norm for an over-constrained system, or to minimize the norm for an under-constrained system. While this makes sense for vectors in which elements are all defined in the same Euclidean coordinate system, it does not make sense when the elements (parameters) are not defined in a common coordinate frame and have incommensurate units. This is because the choice of units affects the magnitude of “error” that is contributed by each element when the sum of squared errors is minimized. For our case the MP inverse attempts to minimize


where the relative contributions of different terms clearly depend on the arbitrary choices of units for the different elements. This shows that squared-error is not a physically meaningful quantity to minimize when different variables are defined in different units. The loss of physical meaning when minimizing Euclidean inner products involving variables with incommensurate units has been observed to be problematic for robotic systems [4, 6], and essentially the same problem applies in our case when the MP inverse is applied to the Jacobian matrix representing the transformation of variables in incommensurate units. This is why our simulation results were not consistent with respect to changes of units from meters to centimeters. The fact that minimizing squared error has no physical meaning in this context reinforces the conclusion that the calculated controls using the MP inverse should not be trusted.

At this point it is valuable to consider what property of a generalized inverse is needed to avoid the unit dependency problems of the MP inverse. What is needed is a solution that is consistent with respect changes of units in the governing equation , which can be expressed in the form of diagonal matrices and as


where is just but with elements defined in different units and is similarly obtained from . The governing equation can now be expressed as


The solution can be expressed in terms of an unknown generalized inverse as


while the solution in the original units was . In order to satisfy our assumption that changes of units take the form of diagonal transformations, e.g., and , the unknown general inverse must satisfy


which is not satisfied by the MP inverse because it only provides consistency with respect to orthogonal transformations (rotations) but not diagonal transformations. The problem to be solved requires a generalized inverse that ensures consistency for diagonal transformations, and an inverse of this kind has been recently developed [10]. This unit-consistent (UC) inverse can be defined in terms of the MP inverse as follows:


where and are positive diagonal matrices determined from , and is a matrix uniquely determined from so that the magnitude of the product of the nonzero elements of each row and column of is 1. (The basis for this decomposition, and its existence and uniqueness properties, are described in [10].) The UC inverse works by eliminating the effect of diagonal transformations. From the definition it can be seen that


So for any diagonal matrices and applied to change the units,


and since and are also diagonal matrices:


This shows that the UC inverse satisfies the consistency requirement of equation 20.

The problems we have demonstrated with the MP inverse occur because it does not produce control behavior that is invariant with respect to the choice of units. Assuming the motion is given as in units of meters, with the Jacobian matrix terms given in equation 9 and , the transformation of units from meters to centimeters, , can be expressed as a diagonal transformation in the form of:


where is the scale factor of converting from meters to centimeters. What is needed is a generalized inverse that will guarantee invariance with respect to units. The UC inverse has this property, so we should expect that using it in place of the MP inverse will eliminate unit dependencies in the control calculations and avoid the unstable behavior demonstrated in our earlier simulations. This can be tested simply by re-running the simulations with the MP inverse replaced by the UC inverse and checking whether the performance of the control system is invariant with respect to changes of units from meters to centimeters. Said another way, changing the units should produce a control solution defined in the new units but which exhibits the exact same physical behavior obtained using the original units.

As was done in the earlier simulations, controls are calculated with a timestep of over a simulation time of for the case of lengths defined in meters and then with lengths in centimeters according to equation 14. Figure 4 shows that in both cases the control produced using the UC inverse converges rapidly, and the velocity variation is identically smooth in both cases. This corroborates our expectation that the UC inverse produces stable results that are not affected by the arbitrary choice of units used for lengths.

It should be noted that we showed in the earlier simulations that the timestep could be reduced so that controls from the MP inverse also produced acceptable results. The difference is that the UC-inverse solution mathematically guarantees that its results are invariant with respect to the choice of units and therefore does not require timestep changes to avoid unpredictable behaviors. The MP inverse does not have the correct mathematical properties and therefore cannot offer the same guarantee.

Fig. 4: Joint velocity () generated by simulations using both length units(meter / centimeter) and timestep (). The result changes smoothly over the simulation time.

Iv Generalized Consistency Considerations

In the previous section it was shown that the UC inverse satisfied all requirements for maintaining unit consistency in the example system333An additional property that is sometimes required of a generalized inverse that was not demonstrated in our example system is consistency with respect to use of the Kronecker product. Because this property was not among those established in [10], we formally prove it in Appendix A.. In this secton we examine how the UC inverse can also be combined with the MP inverse, and even other generalized inverses (e.g., the Drazin inverse [11, 12, 13] or other similarity-consistent inverse [15]), to construct solutions to inverse problems when there is a mix of variables involving different consistency requirements. In addition to variables that require unit consistency to be preserved, a complex real-world system may also involve variables defined in a Cartesian coordinate frame that require consistency with respect to rotations of that coordinate frame. In other words, the behavior of the control system must be invariant with respect to changes of units for some variables and invariant with respect to rotations for other variables. The UC inverse is applicable in one case while the MP inverse is applicable in the other, but what is needed for such a system is a generalized inverse that will guarantee unit consistency for some variables and rotation consistency for others. If we assume444The ordering of the variables is arbitrary so there is no loss of generality in assuming they are permuted so that the UC variables come first and the rotation-consistent variables come next. that the first variables require unit consistency and the remaining variables require rotation consistency then the transformation matrix to be inverted can be block-partitioned as


It has been shown [10] that the mixed inverse can be obtained from this block-partitioned form as


In the previous examples it has been stressed that the problem solution should not depend on arbitrary system preferences because such dependencies lead to unpredictable results. In the case of a system with a mix of consistency requirements we should expect the correct solution to produce the same behavior if units on the first variables are changed or if the coordinate frame of the remaining variables is rotated. If the mixed inverse above is correct then the controls produced should have this property. We will demonstrate that it does in the following example.

Consider the control of a planetary rover with a robotic arm as displayed in Figure 5, which includes two projected views in directions and . The x-y coordinate is shown as frame . The frame , which is an orthogonal transpose of frame by , will be considered later. The rover is free to move in any direction on planar terrain, and its Cartesian position coordinates in this plane are . The part B can rotate and can also ascend/descend555The part B can be thought of as an extendible arm for taking a soil or rock sample. This rover, called ODIF, is intended to function as part of a team of lightweight bots for large area investigation.. In addition, the arm can elongate within a fixed range, but it cannot rotate in the vertical plane. Thus there are 5 degrees of freedom for the design, denoted as , where is the position coordinates of part B.

Fig. 5: Rover with an extendable arm. The rover is free to move on a plane, body part B can ascend/descend or rotate, and the arm can be extended. The two projected views, and , show the structure of the arm. is a fixed angle. The coordinate frame is a rotation of the original frame by an angle of .

From the geometry relations, the position of tip-point is given as ,


We can use the same approach of the previous example to generate the Jacobian matrix for as

The initial states are set to be and , the constant angle . The target velocity of the tip-point is . Since the system has redundant degrees of freedom and therefore is singular, has to be solved with a general inverse (e.g. MP inverse, UC inverse, or mixed inverse). For the unknown variables in , have incommensurate units, and are defined in a common Euclidean space. Thus the Jacobian matrix can be partitioned as

Now the mixed inverse can be used to find the solution for using equation 31. We initially use meters as the length unit and as the coordinate frame. We will then consider a change of length units from meters to centimeters and a coordinate frame rotation from to by a rotation angle of . Given as the scale factor to convert from meters to centimeters, the governing equation for the centimeter and rotated case, , can be expressed as a diagonal transformation of the meter case, , as

This shows the block of the matrix requiring rotation consistency and the block of variables defined in incommensurate units. We then test the three possible approaches to computing the controls: using the MP inverse alone; using the UC inverse alone; and using the mixed inverse obtained from equation 31. The solutions for from the three approaches are displayed in table II for . The column headings give the unit/coordinate frame in which each test was performed but the results are all given in (converted to) a common coordinate frame for comparison purposes. As can be seen, the mixed inverse is the only approach that produces identical results regardless of coordinate-system changes. For the other approaches, it can be seen that when the controls are evaluated with the MP inverse, it generates the same results when different coordinate frames are used but not when length units are changed. When the UC inverse is solely applied to the entire system it generates invariant solutions when units are changed but not when rotations are applied. Therefore, it can be concluded that solely using either the UC inverse or MP inverse alone will not produce reliable results. Instead, the mixed inverse is required to ensure that the behavior of the system is invariant with respect to defined changes of units and coordinates.

A transient simulation was performed for to further observe the full control process. Figure 6(a) displays variation of for the three approaches. It shows that the angular velocity calculated over time by the MP inverse is not affected by a rotation of the coordinate frame from to but is affected by a change of the length unit from meters to centimeters; and the reverse is true for the UC inverse. By contrast, the angular velocity from the mixed inverse is identical over time in all cases.

In summary, for this system involving variables with different consistency requirements the mixed inverse yields reliable control while the alternatives do not. This demonstrates the necessity of using the appropriate inverse to satisfy all applicable consistency requirements.

Length unit() Length unit() Length unit()
Coordinate() Coordinate() Coordinate()
MP inv
UC inv
Mixed inv
TABLE II: Joint velocities , solved with MP inverse, UC inverse, and Mixed inverse approaches.
Fig. 6: Joint velocity for solved with different generalized inverse approaches. (a) MP inverse. (b) UC inverse. (c) mixed inverse. Only the mixed inverse yields the same results over the for the transformations in all three cases.

V Discussion

In this paper we discussed the fact that different generalized inverses have different properties and that the correct one must be chosen based on the properties needed for the problem at hand. For example, the MP inverse is appropriate if the application requires the behavior of the system to be invariant with respect to rotations of the coordinate frame. We then described that many problems require the behavior of the system to be invariant with respect to the choice of units used for key parameters. For example, the behavior should be the same regardless of whether the parameters are defined in metric units or imperial units as long as all parameters are defined consistently with whatever choice of units is made. We emphasized that the Moore-Penrose (MP) pseudoinverse does not satisfy this requirement.

We demonstrated what can happen if the wrong generalized inverse is used by examining the MP inverse in an example of a robotic arm in which consistency with respect to units is required. In that example it was shown that the MP inverse produced different behaviors depending on the choice of units, and as a result it produced erratic and unpredictable behaviors. This is because the MP inverse provides consistency with respect to rotations rather than changes of units. We described that methods are known to tailor such problems by hand so that unit consistency can be maintained, but those methods are complicated and thus are more likely to introduce accidental errors. We explained that what is truly needed to solve such problems is a generalized inverse that is unit consistent. Such an inverse, called the UC inverse, has recently been developed, and we showed that simply replacing the MP inverse with the UC inverse does in fact eliminate all dependencies on the choice of units. Use of the UC inverse led our example system to exhibit stable behavior in all cases.

Lastly we demonstrated that the UC and MP inverses can be combined to provide consistency for systems defined with a mix of variables, some of which demand unit consistency while others require rotational consistency.

Appendix A UC Inverse and the Kronecker Product

The Kronecker product is often used for the mathematical representation of a complex system in terms of simpler subsystems [14]. In this appendix we show that the UC inverse satisfies the same useful properties as the MP inverse with respect to the Kronecker product.

The Kronecker product is a non-commutative tensor operator, usually denoted as

, which takes an matrix and a matrix and constructs a composition of the two matrices to produce a higher-dimensional matrix. The definition of the Kronecker product of and is

Thus each is a of matrix. For example, assuming is and is (i.e., and ) then

The Kronecker product is important in engineering design because it can be used to elegantly and efficiently represent complex systems as compositions of simpler subsystems. It finds applications in control systems, signal processing, image processing, semidefinite programming, and quantum computing[16, 17, 14, 18, 19, 20, 21]. It is bilinear and associative:


and it satisfies the following with respect to the transpose (and conjugate-transpose) operator


For matrices and for which the products and are valid, the following mixed-product property (so-called because it involves both standard matrix multiplication and the Kronecker product) can be shown to hold:


If matrices and are orthogonal then is also orthogonal:


It is also the case that


and more generally


The last two properties are necessary to construct a matrix to transform from one Kronecker-constructed matrix to another Kronecker-constructed matrix, which is required for performing controls of the kinds of robotic and mechanical systems of interest in this thesis, but we also require the ability to apply generalized matrix inverses in the case of singular matrices. It has been proven that the MP inverse satisfies [22]:


and more generally:


but it has not yet been established that these two results also hold for the UC inverse. They must be proven in order to show that the UC inverse can be used for general control systems represented using Kronecker products. We begin by proving that the Kronecker product of two diagonal matrices is also diagonal. Given diagonal matrices and

the Kronecker product

can be seen to also be diagonal because every nonzero element of is at position (), every nonzero element of is at position (), and every nonzero element of is at position . Or more simply, every diagonal block of is a diagonal matrix and therefore must be a diagonal matrix.

With these basic properties in mind, we have the prerequisites to prove the UC inverse property. First consider the base case,


Recall the decomposition of equation 21 for any matrix , where for notational clarity we now use and instead of and :


and the UC generalized inverse of is defined as:


Applying this to the left hand side of equation 46 gives


Using the fact that gives


and the right-hand side of equation 46 is


where and are diagonal matrices. In order to apply equation 48, the rows and columns of must satisfy the product constraint. Let and be represented as

where for any , , , , the matrix elements of and have the following property:


For every row of

the product of its nonzero elements is (when )


The same holds analogously for every column of , so equation 48 can be applied to equation 54 to obtain


and from equation 52 the following theorem can be concluded:

Theorem 1: .

We now show that Theorem 1 can be used as the base case for a mathematical induction proof of the general case involving matrices . Given


it is required to show