Solving a bounded problem

brendanjmeade commented 10 months ago

Summary: Solve the block model problem with a bounded solver so that TDE slip deficit (or fault coupling) rates can be bounded within a specified range. There are three cases to consider:

Classic linear unbounded solve
- Implementation (complete): np.linalg.inv
- Pros:
  - Fast
  - Good for discovery
  - Easy H-matrix implementation for large problems
- Cons:
  - Tendency to potentially oversmooth
  - Coupling values may significantly exceed one without careful hyperparameter tuning

Slip deficit rate bounded solve

Implementation (incomplete): scipy.optimize.lsq_linear

See: celeri_western_north_america_dense_experimental.ipynb

# Example implementation
lower_bound = np.zeros_like(estimation.state_vector)
upper_bound = np.zeros_like(estimation.state_vector)
lower_bound[index.start_tde_col[0] : index.end_tde_col[0]] = 0
upper_bound[index.start_tde_col[0] : index.end_tde_col[0]] = 30
res = lsq_linear(
estimation.operator * np.sqrt(estimation.weighting_vector[:, None]),
estimation.data_vector * np.sqrt(estimation.weighting_vector),
bounds=(lower_bound, upper_bound),
verbose=1,
)
estimation_bounded = copy.deepcopy(estimation)
estimation_bounded.state_vector = res.x
celeri.post_process_estimation(estimation_bounded, operators, station, index)

Pros:
- Works well to recover both block rotation parameters with no bounds and bounded TDE rates
Cons:
- Slower
- Limits discovery
- Not obvious how to H-matrix compress. KL eigenmodes are a better option.

C. Coupling rate bounded solve

Implementation (not started): Iterative use of scipy.optimize.lsq_linear with changing bounds


# Conceptual implementation
coupling_greater_than_one = True
while coupling_greater_than_one:
# Do slip rate bounded estimation
res = lsq_linear(
    estimation.operator * np.sqrt(estimation.weighting_vector[:, None]),
    estimation.data_vector * np.sqrt(estimation.weighting_vector),
    bounds=(lower_bound, upper_bound),
    verbose=1,
)
if max_slip_deficit_tde_rate > fault_slip_deficit_rate:
    # Decrease upper TDE rate bound closer to close to the estimated fault_slip_deficit_rate
else:
    coupling_greater_than_one = False

estimation_bounded = copy.deepcopy(estimation) estimation_bounded.state_vector = res.x celeri.post_process_estimation(estimation_bounded, operators, station, index)



   - Pros:
      - Solve with coupling bounds without having to *a priori* know long-term fault slip rates.
   - Cons:
      - Slower again
      - Limits discovery
      - Not obvious how to H-matrix compress

brendanjmeade commented 4 months ago

CVXOPT and quadratic programming makes this possible and easy. See: celeri_western_north_america_qp.ipynb as a starting point. More do come.

brendanjmeade commented 4 months ago

Proposal for adding slip rate bounds and/or sign constraints

These are nearly identical in practice, so we'll treat sign constraints as bounded between minimum and maximum values.

Specifying bounds in a segment file

Add nine new columns to segment file: ss_rate_bound_flag: 0 or 1. If 0 there is not a bounded strike-slip constraint. If 1 there is a bounded strike-slip constraint. ss_rate_bound_min: minimum strike-slip rate. ss_rate_bound_max: maximum strike-slip rate. ds_rate_bound_flag: 0 or 1. If 0 there is not a bounded dip-slip constraint. If 1 there is a bounded dip-slip constraint. ds_rate_bound_min: minimum dip-slip rate. ds_rate_bound_max: maximum dip-slip rate. ts_rate_bound_flag: 0 or 1. If 0 there is not a bounded tensile-slip constraint. If 1 there is a bounded tensile-slip constraint. ts_rate_bound_min: minimum tensile-slip rate. ts_rate_bound_max: maximum tensile-slip rate.

Need to check on sign convention:

Strike-slip: negative = right-lateral, positive = left-lateral
Dip-slip: negative = CHECK THIS, positive = CHECK THIS
Tensile-slip: negative = CHECK THIS, positive = CHECK THIS

Todo

[ ] Document sign convention (@brendanjmeade)
[ ] Add these to example segment file (@brendanjmeade)
[ ] Add editing capabilities to celeri_ui (@brendanjmeade)
[ ] Develop in QP notebook (@brendanjmeade)
[ ] Clean function call to generate partials (@brendanjmeade)
[ ] Write up linear system for slip rate bounds (@brendanjmeade)
[ ] Write function for generating linear bounds operator (@brendanjmeade)

Thoughts @jploveless ?

jploveless commented 4 months ago

Just want to clarify: are those parameters going into a .command file, or are they new columns of a segment .csv? Seems like the latter, but maybe there are some global settings that need to be applied in a .command?

brendanjmeade commented 4 months ago

Great question. My current plan is new columns of the segment .csv file.

I hadn't thought about global overrides but I can imagine something like command.ignore_all_slip_rate_constraints. Do you have percolating ideas? Happy to consider other options.

jploveless commented 4 months ago

No, that sounds good! I was just clarifying because the last comment said "specifying bounds in a command file".

I'm not opposed to global overrides, but it definitely sounds like something I'd forget to turn off and then scratch my head for a while!

brendanjmeade commented 4 months ago

I totally agree that I'd forget any command overrides as well. I fixed the comment that had the error in it. Thanks!

brendanjmeade / celeri