microsoft/qdk
Publicmirrored from https://github.com/microsoft/qdkAvailable
samples/chemistry/SPSA/src/SPSA.qs
165lines · modecode
| 1 | // Copyright (c) Microsoft Corporation. |
| 2 | // Licensed under the MIT License. |
| 3 | |
| 4 | export |
| 5 | SpsaOptions, |
| 6 | DefaultSpsaOptions, |
| 7 | FindMinimumWithSpsa; |
| 8 | |
| 9 | import Std.Random.DrawRandomInt; |
| 10 | import Std.Arrays.Zipped; |
| 11 | import Std.Arrays.DrawMany; |
| 12 | import Std.Arrays.Mapped; |
| 13 | import Std.Convert.IntAsDouble; |
| 14 | |
| 15 | |
| 16 | /// # Summary |
| 17 | /// Options for use with optimizing objectives via the simultaneous |
| 18 | /// perturbative stochastic approximation (SPSA) algorithm. |
| 19 | /// |
| 20 | /// # Named Items |
| 21 | /// ## StepScale |
| 22 | /// The coefficient by which steps along gradient vectors should be scaled. |
| 23 | /// ## StepPower |
| 24 | /// The power to which the iteration number should be raised when computing |
| 25 | /// how far to step along the gradient vector. |
| 26 | /// ## StepOffset |
| 27 | /// A number to be added to the number of iterations when computing |
| 28 | /// how far to step along the gradient vector. |
| 29 | /// ## SearchScale |
| 30 | /// The coefficient by which searches should be scaled when estimating |
| 31 | /// gradient vectors. |
| 32 | /// ## SearchPower |
| 33 | /// The power to which the iteration number should be raised when computing |
| 34 | /// how far to search in order to estimate gradient vectors. |
| 35 | /// ## NIterations |
| 36 | /// The number of iterations of SPSA to run before stopping. |
| 37 | /// ## MaximumSetback |
| 38 | /// Whether the maximum setback rule is enabled (requiring an additional |
| 39 | /// objective evaluation at each iteration), and if so, the maximum |
| 40 | /// allowed increase in objective values at each iteration. |
| 41 | struct SpsaOptions { |
| 42 | StepScale : Double, |
| 43 | StepPower : Double, |
| 44 | StepOffset : Int, |
| 45 | SearchScale : Double, |
| 46 | SearchPower : Double, |
| 47 | NIterations : Int, |
| 48 | MaximumSetback : (Bool, Double), |
| 49 | } |
| 50 | |
| 51 | /// # Summary |
| 52 | /// Returns a default set of options for use with SPSA optimization. |
| 53 | function DefaultSpsaOptions() : SpsaOptions { |
| 54 | new SpsaOptions { |
| 55 | SearchScale = 0.1, |
| 56 | SearchPower = 0.101, |
| 57 | StepScale = 1.0, |
| 58 | StepPower = 0.602, |
| 59 | StepOffset = 0, |
| 60 | MaximumSetback = (false, 0.1), |
| 61 | NIterations = 30, |
| 62 | } |
| 63 | } |
| 64 | |
| 65 | /// # Summary |
| 66 | /// Given an operation that evaluates an objective at a given point, |
| 67 | /// attempts to find the minimum value of the objective by using the |
| 68 | /// simulntaneous perturbative stochastic approximation (SPSA). |
| 69 | /// |
| 70 | /// # Input |
| 71 | /// ## oracle |
| 72 | /// An operation that evaluates the objective function at a given point. |
| 73 | /// ## startingPoint |
| 74 | /// An initial guess to be used in optimizing the objective function |
| 75 | /// provided. |
| 76 | /// ## options |
| 77 | /// Options used to control the optimization algorithm. |
| 78 | /// |
| 79 | /// # Output |
| 80 | /// The coordinates and final objective value found by the SPSA algorithm. |
| 81 | operation FindMinimumWithSpsa(oracle : (Double[] => Double), startingPoint : Double[], options : SpsaOptions) : (Double[], Double) { |
| 82 | let nParameters = Length(startingPoint); |
| 83 | // The SPSA algorithm relies on projecting gradients onto random vectors |
| 84 | // where each element is either +1 or −1. We can implement that in Q# |
| 85 | // by choosing an element out of [-1.0, +1.0] uniformly at random. |
| 86 | let drawDelta = () => [-1.0, 1.0][DrawRandomInt(0, 1)]; |
| 87 | |
| 88 | mutable currentPoint = startingPoint; |
| 89 | |
| 90 | // Depending on what options are enabled, we may reject certain |
| 91 | // updates, so we keep a counter as to how many iterations have been |
| 92 | // accepted. |
| 93 | mutable nAcceptedUpdates = 0; |
| 94 | mutable lastObjective = 0.0; |
| 95 | |
| 96 | // The SPSA algorithm proceeds by estimating the gradient of the |
| 97 | // objective, projected onto a random vector Δ of ±1 elements. At each |
| 98 | // iteration, the step size used to evaluate the gradient and the |
| 99 | // step taken along the estimated gradient decay to zero, |
| 100 | // such that the algorithm converges to a local optimum by follow |
| 101 | // a directed random walk that is biased by gradients of the objective. |
| 102 | for idxStep in 1..options.NIterations { |
| 103 | Message($"Iteration {idxStep}:"); |
| 104 | |
| 105 | // Following this strategy, we'll start by using the options |
| 106 | // passed into this operation to set αₖ, the amount that we look |
| 107 | // along Δ when using the midpoint formula to evaluate the gradient |
| 108 | // of the objective function 𝑜, and βₖ, the amount that we step |
| 109 | // along the gradient to find the next evaluation point. |
| 110 | let searchSize = options.SearchScale / IntAsDouble(1 + nAcceptedUpdates)^options.SearchPower; |
| 111 | let stepSize = options.StepScale / IntAsDouble(1 + nAcceptedUpdates + options.StepOffset)^options.StepPower; |
| 112 | |
| 113 | // We next draw Δ itself, then use it to find 𝑥ₖ + αₖ Δ and |
| 114 | // 𝑥ₖ − αₖ Δ. |
| 115 | let delta = DrawMany(drawDelta, nParameters, ()); |
| 116 | let search = Mapped(d -> searchSize * d, delta); |
| 117 | let fwd = Mapped((a, b) -> a + b, Zipped(currentPoint, search)); |
| 118 | let bwd = Mapped((a, b) -> a + b, Zipped(currentPoint, Mapped(d -> -d, search))); |
| 119 | |
| 120 | // We then evaluate 𝑜 at each of these two points to find the |
| 121 | // negative gradient 𝑔ₖ = 𝑜(𝑥ₖ − αₖ Δ) − 𝑜(𝑥ₖ + αₖ Δ). |
| 122 | let valueAtForward = oracle(fwd); |
| 123 | let valueAtBackward = oracle(bwd); |
| 124 | let negGradient = (oracle(bwd) - oracle(fwd)) / (2.0 * searchSize); |
| 125 | Message($" obj({fwd}) = {valueAtForward}"); |
| 126 | Message($" obj({bwd}) = {valueAtBackward}"); |
| 127 | |
| 128 | // We can step along 𝑔ₖ to find 𝑥ₖ₊₁. Depending on whether options |
| 129 | // such as the maximum setback rule are enabled, we may reject |
| 130 | // the update. Either way, we report out to the caller at this |
| 131 | // point. |
| 132 | let step = Mapped(d -> negGradient * stepSize * d, delta); |
| 133 | let proposal = Mapped((a, b) -> a + b, Zipped(step, currentPoint)); |
| 134 | if Fst(options.MaximumSetback) { |
| 135 | // Is this our first update? If so, accept and set the |
| 136 | // lastObjective. |
| 137 | if nAcceptedUpdates == 0 { |
| 138 | Message($" First update; accepting."); |
| 139 | lastObjective = oracle(proposal); |
| 140 | nAcceptedUpdates += 1; |
| 141 | currentPoint = proposal; |
| 142 | } else { |
| 143 | // How much did our objective get worse (increase) by? |
| 144 | let thisObjective = oracle(proposal); |
| 145 | if thisObjective - lastObjective <= Snd(options.MaximumSetback) { |
| 146 | Message($" Proposed update gave objective of {thisObjective}, which is within maximum allowable setback of previous objective {lastObjective}; accepting."); |
| 147 | // Within the limit, so we're good. |
| 148 | lastObjective = thisObjective; |
| 149 | nAcceptedUpdates += 1; |
| 150 | currentPoint = proposal; |
| 151 | } else { |
| 152 | Message($" Proposed update gave objective of {thisObjective}, which exceeds maximum allowable setback from previous objective {lastObjective}; rejecting."); |
| 153 | } |
| 154 | } |
| 155 | } else { |
| 156 | // No maximum setback rule, so always accept the proposed |
| 157 | // update. |
| 158 | nAcceptedUpdates += 1; |
| 159 | currentPoint = proposal; |
| 160 | } |
| 161 | |
| 162 | } |
| 163 | |
| 164 | return (currentPoint, oracle(currentPoint)); |
| 165 | } |
| 166 | |