Update factors of a matrix

New to Julia
1

I now think I understand how Julia modifies its arguments (by changing what the arguments contain) and I have been doing so using .=, as in, e.g.:

function foo!(x)
  x .= x + 1
end

But what if the argument is more complicated, i.e., if I want to modify the factors of a matrix: Say that I have a variable J which contains the factors of a 3x3 matrix A:

J = factorize(A)

And then I would like to write a function foo! that replaces what J contains with some new factors, that is, something like:

function foo!(J)
   # Create a 3x3 matrix B
   J .= factorize(B)
end

This does not work because J is not a simple array. So how does one write a function foo! that modifies the content of J (so that it contains the output of factorize(B))?

2

Update: I have tried

J.factors .= factorize(B).factors

which seems to work… I am unsure it is the correct way though…

3

Factorizations are implemented as structs, eg

julia> dump(factorize(randn(3,3)))
Base.LinAlg.LU{Float64,Array{Float64,2}}
  factors: Array{Float64}((3, 3)) [-1.29583 1.08384 0.832042; -0.882121 2.24427 0.94113; -0.457674 -0.310311 2.52631]
  ipiv: Array{Int64}((3,)) [1, 2, 3]
  info: Int64 0

I am not sure that updating them in place is worthwhile; except for a few special applications, the cost of the factorization will dominate the allocation. But if you insist, you should probably write specific methods for each factorization type you deal with.

4

Oh I see… Thank you for your reply. However, this does not help with my specific problem. Maybe I should detail it a bit more. I have to (on some condition) update and factorize a large n\times n sparse matrix J, which represents the Jacobian J(x) = \frac{\partial f}{\partial x}(x) of a function f:\mathbb{R}^n\rightarrow\mathbb{R}^n of a variable x \in \mathbb{R}^n. I have already coded a function for the jacobian:

function Jacobian(x,...)
  # some calculations
  return J
end

Then I have two functions, foo and bar, that require the updated Jacobian J. In MATLAB I would have used the (frowned upon) global variable declaration, so foo and bar would look like

function foo(x,...) # or function bar(x,...)
  global J
  if J_needs_update(x,...)
    J = factorize(Jacobian(x))
  end
  # some calculations that use the factorized J
  return ...
end

This would ensure that an updated J is always available to foo and bar. If foo updates J, then bar does not need to update it, and vice versa. (And the costly factorization is only done when necessary.) I understand that global variables is generally not the correct way to do things. However if I use J as an argument to foo and bar which do not modify its value, then foo and bar do not communicate the updated (and factorized) J anymore. How does one do what I want to do without the use of a global variable? Ideally, I would just have foo! and bar! modify J, as in, e.g.,

function foo!(x,J,...) # or function bar!(x,J,...)
  if J_needs_update(x,...)
    J = factorize(Jacobian(x)) # but modify the argument J 
  end
  # some calculations that use the factorized J
  return ...
end

(Please let me know if something is unclear)

On a side note, I thought it would be pretty neat if appending a ! at the end of a function was not just a convention, but added the property that any argument that is modified within the function is modified outside. What I mean is for example that

function add_one!(x)
  x = x + 1
end

actually modifies a if add_one!(a) is called, just because add_one! has a ! appended, even though I did not use x .= x + 1 inside. But maybe this is a stupid remark, I’m not sure…

5

This is a perfect use case for defining your own (mutable) struct to hold the factorization. For example, you can simply do:

mutable struct MyFactorization
  J
  needs_update::Bool
end

This struct is mutable, so you can pass it into a function and do:

function foo(fact::MyFactorization, x)
  if fact.needs_update
    fact.J = factorize(Jacobian(x))
    fact.needs_update = false
  end
end

Any mutation of fact inside the function will be visible to anywhere else you have that same fact object. So if you do:

fact = MyFactorization(...whatever...)
foo(fact, x)
@show fact.J

then you’ll see that fact.J will have been updated. This avoids any global variables and means you could even have multiple factorization objects, should you need that in the future.

For better performance, you can annotate the type of J in your struct:

mutable struct MyFactorization
  J::LinAlg.LU{Float64,Array{Float64,2}}  # I'm guessing here, your type may be different
  needs_update::Bool
end

Or for better performance and more flexibility, you can make that type a parameter, and Julia will figure it out for you:

mutable struct MyFactorization{T}
  J::T
  needs_update::Bool
end

Even better, once you’ve done this, you can make some convenient functions, like:

function factors(f::MyFactorization, x)
  if f.needs_update
    f.J = factorize(Jacobian(x))
    f.needs_update = false
  end
  return f.J
end

and now you can simply call factors(f, x) everywhere and it will just do the right thing, only updating the factorization as necessary.

On a side note, I thought it would be pretty neat if appending a ! at the end of a function was not just a convention, but added the property that any argument that is modified within the function is modified outside.

That would be interesting, but it doesn’t really fit into the way variables are passed in Julia (it’s more like what’s done with references in C++), and it would spectacularly break a lot of existing code, so it’s not likely to happen any time soon :wink:

1 Like