PDLPorters / pdl-linearalgebra

https://p3rl.org/PDL::LinearAlgebra
1 stars 7 forks source link

The matrix multiplication 'x' drops the imaginary component, but matmult works. #18

Closed KJ7LNW closed 2 years ago

KJ7LNW commented 2 years ago

here is the setup with trivial values using the latest PDL from git (8c46144b9f6a2788673a78f7697eeb7f764c8def):

pdl> use PDL::LinearAlgebra

pdl>  $S = sequence(2,2) + sequence(2,2)*i

pdl> p $sqrt_y = stretcher(zeroes(2)+sqrt(50))->minv

[
 [0.14142136          0]
 [        -0 0.14142136]
]

pdl> p $id = identity(2)

[
 [1 0]
 [0 1]
]

pdl> p $id - $S

[
 [    1  -1-i]
 [-2-2i -2-3i]
]

Here is the inconsistency after the setup above; note that the x operator has no imaginary part, but matmult does:

pdl> p $sqrt_y x ($id-$S)

[
 [ 0.14142136 -0.14142136]
 [-0.28284271 -0.28284271]
]

pdl> p matmult($sqrt_y, ($id-$S))

[
 [                     0.14142135623731   -0.14142135623731-0.14142135623731i]
 [-0.282842712474619-0.282842712474619i -0.282842712474619-0.424264068711929i]
]

Here are the real-life numbers from a 10pF capacitor:

setup:

pdl> use PDL::LinearAlgebra

pdl> p $S = pdl [ [0.7163927-0.4504119*i, 0.2836073+0.4504119*i ], [0.2836073+0.4504119*i, 0.7163927-0.4504119*i ] ] 
[
 [0.7163927-0.4504119i 0.2836073+0.4504119i]
 [0.2836073+0.4504119i 0.7163927-0.4504119i]
]

pdl> 
pdl> p $sqrt_y = stretcher(zeroes(2)+sqrt(50))->minv

[
 [0.14142136          0]
 [        -0 0.14142136]
]

pdl> p $id = identity(2)

[
 [1 0]
 [0 1]
]

pdl> p $id-$S
[
 [ 0.2836073+0.4504119i -0.2836073-0.4504119i]
 [-0.2836073-0.4504119i  0.2836073+0.4504119i]
]

and the problem with missing i values.

pdl> p $sqrt_y x ($id-$S)

[
 [ 0.040108129 -0.040108129]
 [-0.040108129  0.040108129]
]

pdl> p matmult($sqrt_y, ($id-$S))

[
 [ 0.0401081290048015+0.0636978617634234i -0.0401081290048015-0.0636978617634234i]
 [-0.0401081290048015-0.0636978617634234i  0.0401081290048015+0.0636978617634234i]
]

The matmult numbers provide the same thing I get in Math::Matrix::Complex so I think the numbers are right, its just an oddity with 'x'.

KJ7LNW commented 2 years ago

More info:

This is a PDL::LinearAlgebra problem because it works natively:

pdl> p $S = pdl [ [0.7163927-0.4504119*i, 0.2836073+0.4504119*i ], [0.2836073+0.4504119*i, 0.7163927-0.4504119*i ] ] 
[
 [0.7163927-0.4504119i 0.2836073+0.4504119i]
 [0.2836073+0.4504119i 0.7163927-0.4504119i]
]

pdl> 
pdl> p $sqrt_y = stretcher(zeroes(2)+sqrt(50))->inv

[
 [0.14142136          0]
 [         0 0.14142136]
]

pdl> p $id = identity(2)

[
 [1 0]
 [0 1]
]

pdl> 
pdl> p $id-$S

[
 [ 0.2836073+0.4504119i -0.2836073-0.4504119i]
 [-0.2836073-0.4504119i  0.2836073+0.4504119i]
]

pdl> p $sqrt_y x ($id-$S)

[
 [ 0.0401081290048015+0.0636978617634234i -0.0401081290048015-0.0636978617634234i]
 [-0.0401081290048015-0.0636978617634234i  0.0401081290048015+0.0636978617634234i]
]
mohawk2 commented 2 years ago

This is explicitly a PDL::LinearAlgebra problem and should be moved to there (am doing so now).

mohawk2 commented 2 years ago

I believe https://github.com/PDLPorters/pdl-linearalgebra/blob/a162b361f076f5bd93a363a784af0df92ef90ef0/Real/real.pd#L91-L97 is only checking the left arg for being complex, and it should instead check both; the overload's conditional operator condition needs to be !(grep ref($_) && !$_->type->real, @_) instead.

mohawk2 commented 2 years ago

Now released as 0.35.