"For some reason, inv(x'x)x'y != inv(x'x)*x'*y"
x'x is lowered to Ac_mul_B(x, x)
So the left hand side can be rewritten as:
inv(Ac_mul_B(x, x)) * Ac_mul_B(x, y) == inv(x'x)x'y
The right hand side:
inv(Ac_mul_B(x, x)) * transpose(x) * y == inv(x'x)*x'*y
So the order in which the multiplications are evaluated are different and
the result will be subject to normal floating point inaccuracies.
On Wednesday, May 25, 2016 at 1:52:59 PM UTC+2, Zanna Iscenko wrote:
>
> Dear all
>
> I am quite new to Julia, and I am experiencing very odd failures of
> associative matrix multiplication results depending on the permutation of
> the inv(), multiplication signs and whether or not an intermediate
> multiplication result is defined as a separate variable or not.
>
> I'm going to provide a replication with a small toy example later, but
> here is an outline of the actual regression results where this first came
> up (happy to provide data for the underlying matrices if needed)
>
> x_dum : 104x34 Array{Real, 2}
> delta0: Vector, float64, 104
>
> julia> Px_dum = inv(x_dum'*x_dum)*x_dum'
> julia> theta10 = Px_dum*delta0
> (Vector, Real, 34)
> julia> theta11 = inv(x_dum'*x_dum)*(x_dum'*delta0)
> (Vector, Real, 34)
>
> julia>theta10[1:8]
> 8-element Array{Real,1}:
> 0.863612
> 0.315274
> -0.18194
> -0.99634
> -3.16722
> -2.40826
> -1.78346
> -4.2834
>
> julia>theta11[1:8]
> 8-element Array{Real,1}:
> 9.0
> 0.315274
> -0.18194
> -0.99634
> -3.25
> -7.6875
> -5.9375
> -7.4375
>
> I'm genuinely at my wits' end from trying to determine whether this is
> arising from me overlooking some basic (if counter-intuitive) Julia syntax
> rules or this is a genuine quite profound bug.
> All this has been observed while running version 0.4.5 on Windows 64bit,
> and is reproducible in REPL and scripts pre-written/run through julia IDE
> in Atom.
>
> Now, for a toy example to demonstrate when I encounter this result more
> generally.
>
> julia> x=[1 2; 3 4]
> julia> y=[3 4]'
>
> Let's try a basic regression model in different formulations:
>
> julia> inv(x'x)x'y ==inv(x'x)*(x'y)
> true
> julia> inv(x'x)*x'y ==inv(x'x)*(x'y)
> true
> julia> inv(x'x)*x'y ==inv(x'x)*(x'*y)
> true
> julia> inv(x'x)*x'*y ==inv(x'x)*(x'*y)
> false
> ??? Why would this be happening? (the same occurs with inv(x'*x) in the
> first term)
>
> For some reason, inv(x'x)x'y != inv(x'x)*x'*y and inv(x'x)*(x'*y) !=
> inv(x'x)*x'*y This is pretty insidious because it's quite common to use
> parentheses to keep track of variables and/or * to separate them out.
>
> I've tried to trace in what circumstances the issue appears:
>
> 1. Is it the inv() call that is messing with things? No.
> julia> z=inv(x'x)
> julia> z*x'*y == z*x'y
> false
>
> (This does not occur if z is some other matrix)
>
> 2. Can the issue be avoided through using \ rather than inv? No.
> julia> z = x'x
> julia> inv(z)*x'*y == inv(z)*(x'*y)
> false
> julia> z\x'*y == z\(x'*y)
> false
>
> 3. Am I messing up because I combine transposition with * (x'* sequence)?
> No.
>
> julia> Px = inv(x'x)x'
> julia> Px*y == inv(x'x)x'y
> false (whoa!)
>
> BUT
> julia> Px*y ==inv(x'x)x'*y== inv(x'x)*x'*y
> true
>
> 4. Does it happen when y is a matrix instead of a vector? Yes. (Unless y
> is an identity matrix)
> julia> y = [1 5; 8 5]
> julia> inv(x'x)x'y == inv(x'x)x'*y
> false
>
> So, the issue seems to be about the multiplication sign before y even
> though it makes no difference in other contexts.
>
> In small examples the difference might be just floating point noise (e.g.
> try
> inv(x'x)x'y and inv(x'x)*x'*y for x=[5 3; 9 10] and y=[1000; 1000]).
> But as shown above, for real-life cases the divergence can become very
> large.
>
> Would be incredibly grateful for any insight on why this might be
> happening and whether it's actually a bug that needs to be reported.
> Apologies if this has come up before: I tried searching the message board
> and julia issues list on GitHub and didn't find anything relevant.
>
> All best
>
> Zanna
>
>
>
