Ax=B solving for non-square A-CodePudding

I am having trouble solving the following problem. `


# This generates the target Matrix

A =  matrix(c(0.51,0.42,   0,   0,   0,   0,   0,   0,   0,   0,
              0   ,   0,0.35,0.49,0.43,0   ,   0,   0,   0,   0,
              0   ,   0,   0,   0,   0,0.43,0.31,0.97,   0,   0,
              0   ,   0,   0,   0,   0,   0,   0,   0,0.70,0.42,
              1   ,   0,   1,   0,   0,   0,   0,   0,   0,   0,
              0   ,   1,   0,   1,   0,   1,   0,   0,   0,   0,
              0   ,   0,   0,   0,   1,   0,   1,   0,   1,   0,
              0   ,   0,   0,   0,   0,   0,   0,   1,   0,   1),
            nrow = 8, byrow = T)


# x is equal to the coefficients vector 

x = c(50,36,60,85,14,22,84,92,34,74)


# y is equal to the product Ax = y

y = A %*% X

So let's say I did not know A but I do know x and y. How can I find the original Matrix A (or at very least a Matrix which would satisfy Ax=y) from only the x and y vectors?

I really appreciate your time, thank you in advance.

I have tried to follow the equation

A = BX^-1

However, this lead nowhere and the issue I think is around the non-square nature of the A matrix.

CodePudding user response：

The problem is that there is more than one non-square matrix A that satisfies A %*% x == y when x and y are known, so there is no way to take x and y and get your specific matrix back.

However, it is possible to take x and y and return a matrix that satisfies the formula. It requires finding the generalized inverse of x:

library(MASS)

z <- ginv(x)

z
#>             [,1]         [,2]        [,3]        [,4]
#> [1,] 0.001335007 0.0009612047 0.001602008 0.002269511
#>              [,5]         [,6]        [,7]        [,8]
#> [1,] 0.0003738018 0.0005874029 0.002242811 0.002456412
#>              [,9]      [,10]
#> [1,] 0.0009078044 0.00197581

If we multiply y by z then we will get a matrix A that satisfies the condition A %*% x = y

A <- y %*% z

A
#>            [,1]       [,2]       [,3]       [,4]       [,5]
#> [1,] 0.05422797 0.03904414 0.06507356 0.09218754 0.01518383
#> [2,] 0.09167490 0.06600593 0.11000988 0.15584733 0.02566897
#> [3,] 0.16652872 0.11990068 0.19983446 0.28309882 0.04662804
#> [4,] 0.07326516 0.05275091 0.08791819 0.12455077 0.02051424
#> [5,] 0.14685072 0.10573252 0.17622086 0.24964622 0.04111820
#> [6,] 0.19090594 0.13745227 0.22908712 0.32454009 0.05345366
#> [7,] 0.17622086 0.12687902 0.21146504 0.29957547 0.04934184
#> [8,] 0.22161109 0.15955998 0.26593330 0.37673885 0.06205110
#>            [,6]       [,7]       [,8]       [,9]      [,10]
#> [1,] 0.02386030 0.09110298 0.09977946 0.03687502 0.08025739
#> [2,] 0.04033696 0.15401383 0.16868181 0.06233893 0.13567885
#> [3,] 0.07327264 0.27976824 0.30641284 0.11323953 0.24646250
#> [4,] 0.03223667 0.12308547 0.13480789 0.04982031 0.10843244
#> [5,] 0.06461432 0.24670921 0.27020532 0.09985849 0.21733906
#> [6,] 0.08399861 0.32072197 0.35126692 0.12981604 0.28254078
#> [7,] 0.07753718 0.29605105 0.32424639 0.11983019 0.26080688
#> [8,] 0.09750888 0.37230662 0.40776440 0.15069554 0.32798441

Now, testing we have:

A %*% x
#>        [,1]
#> [1,]  40.62
#> [2,]  68.67
#> [3,] 124.74
#> [4,]  54.88
#> [5,] 110.00
#> [6,] 143.00
#> [7,] 132.00
#> [8,] 166.00

This is (within the limits of floating point arithmetic) the same as y:

all(abs((A %*% x) - y) < 1e-12)
#> [1] TRUE

CodePudding user response：

There is no unique solution for A (and no solution at all if x is 0 and y is non zero) but one possible A, call it AA, is the following provided x has no zero elements.

AA <- y %*% (1 / (length(x) * x))

all.equal(AA %*% x, y)
## [1] TRUE

If x can have 0 elements but is not entirely 0 then use this more general formula

AA <- y %*% (1 / (sum(x != 0) * replace(x, x == 0, 1)))
all.equal(AA %*% x, y)
## [1] TRUE