A good match to an assay does not necessarily identify a unique mineral mixture. This worked synthetic example shows how to distinguish what the chemistry supports, what depends on assumptions, and which additional measurement could improve the decision.
Read the interpretation first; expand the code when you want to examine the implementation.
Why bother converting elements to minerals?
What interests me about element-to-mineral conversion (EMC) is fairly simple: we usually have far more assay data than mineralogical measurements. The question is how much mineralogical information we can recover from it, and what we have to assume along the way.
Sometimes the assay gets us a useful distance. Sometimes it gives us a mineral total but cannot tell us how that total is split between phases. That distinction matters before we start using the output to explain metallurgical behaviour.
An assay measures chemistry. EMC estimates mineral proportions using that chemistry, a candidate mineral list and the composition assigned to each mineral. It does not directly measure those proportions, and bulk chemistry alone does not recover texture or liberation. The mineral hosting an element may matter to processing, but a mineral estimate is still not a recovery model.
I think the useful question is not whether an algorithm can return a mineral composition. It is whether the data support the mineral distinctions we need for the decision.
This post works through one invented mixture to make that question concrete. The calculations are original teaching examples, not reproductions of the datasets or complete algorithms in the papers discussed below. Everything uses mass fractions internally; 0.25 means 25 wt%, or 25 g in a 100 g sample.
What the papers contribute
EMC is not a new idea. The useful progression is from solving a chemical balance to asking what makes that balance identifiable, geologically plausible and useful at a particular deposit.
Johnson, Chu and Hussey (1985): the coefficients need checking too. Their simultaneous-equation approach combined chemical and physical properties to estimate clay-mineral proportions. Candidate minerals, plausible component properties and residual checks were all part of the workflow. I take a practical point from this: the mineral composition matrix is an assumption to investigate, not just a table to feed into a solver.
Whiten (2008): an exact fit can have several answers. Singular-value decomposition provides a way to examine which mineral distinctions the chemical system can resolve. This is the starting point for the example below. A pseudoinverse can return one solution while a whole family of other mineral mixtures gives exactly the same assay.
Berry, Hunt and McKnight (2011), then Lund, Lamberg and Lindberg (2013): bring in measurements that address the ambiguity. Berry and colleagues discuss using sparse QXRD to calibrate an EMC calculation, and combining chemistry and QXRD through error-weighted fitting. They also emphasise water and carbon dioxide measurements. At Malmberget, Lund and colleagues used chemical analyses together with mineral chemistry and magnetic measurements, with mineralogical comparisons. That is an important operational example: useful conversion came from a characterised mineral system and supporting measurements, not a more elaborate solver alone.
Parian and colleagues (2015), and Mena Silva and colleagues (2018): calibration is system-specific. Parian and colleagues combined EMC and quantitative XRD for estimating mineral grades. Mena Silva and colleagues compared least-squares and regression-based conversions at the Nabbaren nepheline-syenite deposit; their regression approach performed better in that study. Neither result establishes a universally transferable conversion matrix. The labelled mineral system, analytical methods and calibration data remain part of the model.
Rodrigues and colleagues (2023): make the uncertainty explicit. Their Pedras conference abstract describes Bayesian approximation of modal mineralogy from assays. It supports that broad direction, but an abstract is not enough to reproduce its priors, likelihood or validation. The Bayesian example here is therefore my own deliberately small model, not an implementation of Pedras.
These papers also help separate two ideas that are often bundled together. A normative calculation allocates chemistry according to a prescribed scheme. An EMC model can instead aim to estimate the actual mineral assemblage using independently identified phases. Both involve assumptions, but an estimate of modal mineralogy is not automatically a CIPW-style normative mineralogy (Lund, Lamberg and Lindberg, 2013).
The parts I implement below are the forward mass balance, null-space analysis, constrained fitting, feasible bounds and a conditional Bayesian update. I do not reproduce the papers’ deposit-specific calibrations or claim to implement HSC’s internal algorithms.
One mixture, first forwards
Imagine 60 g quartz, 25 g anatase and 15 g rutile in a 100 g sample. Anatase and rutile are polymorphs: different crystal structures with the same ideal formula, TiO₂. In this simplified oxide basis, quartz contributes SiO₂, while both titanium phases contribute TiO₂.
Let \(m\) contain the mineral mass fractions and \(a\) contain the measured oxide mass fractions. The response matrix \(C\) has assay components in its rows and minerals in its columns:
Entry \(C_{ij}\) is the amount of assay component \(i\) contributed by one unit mass of mineral \(j\). The multiplication is just a weighted sum, \(a_i=\sum_j C_{ij}m_j\). Here it gives
where \(n_{ij}\) is the number of atoms of element \(i\) in mineral \(j\) and \(M_i\) is its atomic molar mass. Convert an elemental response to an oxide-equivalent response using \(M_{\mathrm{oxide}}/(\nu M_{\mathrm{element}})\), where \(\nu\) counts the relevant element in that oxide.
The sample and matrix must use the same reporting basis. Converting FeO to Fe₂O₃-equivalent changes the reporting convention; it does not measure iron oxidation state. Mineral chemistry can also differ from the ideal formula.
Whole-sample mineral closure, \(\sum_jm_j=1\), is not the same as normalising the reported assays. With incomplete assay coverage, water, carbon dioxide or oxide-equivalent reporting, the columns of \(C\) need not sum to one. If retained assay channels have been normalised, their forward model is instead
provided \(C\) represents those same channels before normalisation. Check the dry basis, volatiles and normalisation before fitting anything. Closing an incomplete mineral list to 100% is also an assumption, not a correction for missing phases.
Why the reverse calculation has several answers
Starting with the assay, we want to recover \(m\). There is no general operation called “divide by \(C\)”. An ordinary inverse exists only for a square, full-rank matrix. This one has two rows and three columns.
More importantly, the anatase and rutile columns are identical. Every mixture
This is a null direction. The matrix has rank two, leaving one unconstrained direction among its three mineral fractions. The pseudoinverse gives \((0.60,0.20,0.20)^T\): the minimum-Euclidean-norm solution among the least-squares solutions. It does not discover that the sample contained equal amounts of the two polymorphs (Whiten, 2008; NumPy Developers, n.d.).
Every point on these lines gives the same 60/40 assay. The total titanium phase is identified; its split is not.
Closure needs a qualification here. In this example the closure row is already the sum of the two assay rows. Adding \(\mathbf1^Tm=1\) does not increase the rank. In another mineral system it may add an independent constraint. Non-negativity can also make a particular boundary solution unique. So “more minerals than assays” is a warning, not a complete identifiability test.
For the equality constraints, the general question is whether a non-zero change \(v\) can satisfy both \(Cv=0\) and \(\mathbf1^Tv=0\). Full column rank of the augmented matrix \([C;\mathbf1^T]\) rules out such changes. When it does not have full column rank, inspect which changes remain feasible under non-negativity and any other justified constraints.
For this mixture the answer is straightforward: quartz and total titanium phase are identifiable, but the polymorph split is not. Better precision in the same two channels cannot separate exact duplicate columns. Nearly identical columns present a related but different problem: the split can be identifiable in principle yet very sensitive to noise.
What the solvers actually choose
With noisy assays and a fixed response matrix, a common model is
where \(W^TW=\Sigma_a^{-1}\) when the covariance is positive definite. For independent assay errors, \(W\) has entries \(1/\sigma_i\) on its diagonal. Those are measurement-error weights, not a device for forcing a preferred mineral answer. A fully normalised assay can have singular covariance, so blindly inverting its covariance is not appropriate.
The exact toy example below uses equal weights, not a fitted uncertainty model. NNLS adds non-negativity but does not impose closure. The separate constrained fit imposes both. Their mineral splits need not agree, even when their assay fits do (SciPy Community, n.d.).
Linear programming gives a different answer to a different question. For each mineral, minimise and maximise its fraction subject to the exact assay balance, non-negativity and closure. That finds the feasible range rather than selecting a preferred composition.
View code
m_nnls, _ = nnls(C, a)fit = minimize(lambda m: 0.5* np.sum((C @ m - a) **2), x0=np.full(3, 1/3), jac=lambda m: C.T @ (C @ m - a), method="SLSQP", bounds=[(0., 1.)] *3, constraints={"type": "eq", "fun": lambda m: m.sum() -1.,"jac": lambda m: np.ones(3)}, options={"ftol": 1e-12, "maxiter": 1000},)ifnot fit.success:raiseRuntimeError(fit.message)bounds_lp = []for j inrange(3): endpoints = []for sign in (1., -1.): result = linprog(sign * np.eye(3)[j], A_eq=M, b_eq=np.r_[a, 1.], bounds=(0., None), method="highs")ifnot result.success:raiseRuntimeError(result.message) endpoints.append(result.x[j]) bounds_lp.append(endpoints)bounds_lp = np.asarray(bounds_lp)np.testing.assert_allclose(bounds_lp, [[0.60, 0.60], [0., 0.40], [0., 0.40]], atol=atol)solver_rows = []for label, m in {"Pseudoinverse": m_pinv, "NNLS": m_nnls,"Non-negative + closure": fit.x}.items(): np.testing.assert_allclose(C @ m, a, atol=atol)assert m.min() >=-atol np.testing.assert_allclose(m.sum(), 1., atol=atol) # True for this toy only. solver_rows.append({"Method": label, **dict(zip(names, 100* m)),"Residual norm": np.linalg.norm(C @ m - a)})display(DataFrame(solver_rows).round(6))display(DataFrame({"Mineral": names, "Minimum wt%": 100* bounds_lp[:, 0],"Maximum wt%": 100* bounds_lp[:, 1]}))
Method
Quartz
Anatase
Rutile
Residual norm
0
Pseudoinverse
60.0
20.0
20.0
0.0
1
NNLS
60.0
40.0
0.0
0.0
2
Non-negative + closure
60.0
20.0
20.0
0.0
Mineral
Minimum wt%
Maximum wt%
0
Quartz
60.0
60.0
1
Anatase
0.0
40.0
2
Rutile
0.0
40.0
The result is 60 wt% quartz and 0–40 wt% for either titanium mineral, with their total fixed at 40 wt%. These are feasibility bounds, not confidence or credible intervals. The fractions are coupled: 40 wt% anatase and 40 wt% rutile cannot occur together.
The residuals are all effectively zero. Plotting differences at floating-point scale would not tell us which mineral split is better. Nor should an exact feasible-set calculation be carried over unchanged to noisy data: specify an error model or justified measurement bounds rather than quietly loosening a solver tolerance.
Allocation rules, search and calibration
Staged allocation can be useful when the rules have a mineralogical basis. But it is important to separate the order of the calculation from the rule choosing a mineral. In this toy, allocating SiO₂ before or after TiO₂ makes no difference if the titanium recipient stays fixed. Assigning TiO₂ to anatase rather than rutile changes the answer because the tie-break changed.
View code
allocation_rows = []for ti_phase in (1, 2): finals = []for order in ((0, 1), (1, 0)): m = np.zeros(3) residual = a.copy()for component in order: mineral =0if component ==0else ti_phase amount = residual[component] / C[component, mineral] m[mineral] += amount residual -= amount * C[:, mineral] np.testing.assert_allclose(C @ m, a, atol=atol) finals.append(m) allocation_rows.append({"Order": "Si first"if order[0] ==0else"Ti first","Ti recipient": names[ti_phase],**dict(zip(names, 100* m))}) np.testing.assert_allclose(finals[0], finals[1], atol=atol)display(DataFrame(allocation_rows))
Order
Ti recipient
Quartz
Anatase
Rutile
0
Si first
Anatase
60.0
40.0
0.0
1
Ti first
Anatase
60.0
40.0
0.0
2
Si first
Rutile
60.0
0.0
40.0
3
Ti first
Rutile
60.0
0.0
40.0
A genetic algorithm also searches the model we give it. For this fixed linear response matrix, the squared-error objective with simplex constraints is convex; a genetic search is unnecessary. It may become relevant to a genuinely nonlinear or discrete model, but changing its seed is not a way to quantify posterior uncertainty. A candidate list containing only quartz and anatase can fit our assay exactly while excluding the rutile that we put into the sample.
Calibration is different because it uses additional data. For example, assigning weights of 0.70 and 0.30 to two endpoint solutions gives
Those weights are illustrative choices here. No labels were used to fit them. A real calibration would need the paired data, fitting objective and separate evaluation to be shown. It should not be described as a reproduction of Berry’s or Mena Silva’s methods.
The same distinction applies to a neural model: modal labels can train an assay-to-mineral mapping, and a mass-balance layer can check reconstructed chemistry. The reconstruction check is not an additional mineral measurement. A learned mapping may predict a typical split in its training population without identifying the split in an individual chemically indistinguishable sample.
What would actually separate the two minerals?
An additional measurement needs to change along the unresolved direction. If a new linear observation has response row \(h\), then it distinguishes this direction when
\[
hv\ne0.
\tag{4}\]
For a direct anatase-sensitive response, \(h=(0,1,0)\) and \(hv=1\). Another total-TiO₂ measurement has \(h=(0,1,1)\) and \(hv=0\). One separates the alternatives; the other does not. The noise level then determines how precisely the split can be estimated.
This is the practical link back to the papers: choose mineralogical, mineral-chemical or other diagnostic measurements that address the ambiguity in the actual mineral system. Do not assume that any additional measurement will help equally.
Bayesian inference: uncertainty about which quantity?
Bayesian inference makes the assumptions and uncertainty explicit. It does not make identical chemical signatures different.
A Dirichlet prior is one option for a closed mineral composition:
All \(\alpha_j\) must be positive. The Dirichlet also imposes negative pairwise covariances; that is a restriction of this distribution, not a universal property of compositional data. It has no point mass at an exactly absent mineral. These are modelling choices worth checking, especially for mineral assemblages containing genuine absences (SciPy Community, n.d.).
Condition on the total before discussing the split
Suppose initially \(m\sim\operatorname{Dirichlet}(6,2,2)\). Define the total titanium phase and its anatase share as
\[
S=m_A+m_R,\qquad u=\frac{m_A}{S}.
\]
Under this prior, \(S\) and \(u\) are independent, with \(S\sim\operatorname{Beta}(4,6)\) and \(u\sim\operatorname{Beta}(2,2)\). Our exact toy assay fixes \(S=0.40\) and quartz at 0.60. Conditional on that total,
The factorisation can be checked directly. Set \(m_Q=1-S\), \(m_A=Su\) and \(m_R=S(1-u)\). The absolute Jacobian is \(S\), so the transformed Dirichlet density is proportional to
\[
S^{2+2-1}(1-S)^{6-1}\,u^{2-1}(1-u)^{2-1}.
\]
It separates into a Beta density for \(S\) and another for \(u\).
The \(1/0.40\) factor accounts for the change of scale. This is the conditional prior along our feasible line, not the original marginal prior for anatase, which was \(\operatorname{Beta}(2,8)\).
That distinction matters. The assay has informed the total. It has not informed the share \(u\) within that total. Only in this conditional, fixed-total problem is the assay likelihood constant across \(t\), leaving the conditional prior unchanged. With noisy totals or other prior dependencies, mineral-fraction marginals can still update indirectly.
Add an invented phase-sensitive observation
Now introduce a synthetic observation of anatase on the same mass-fraction scale:
\[
y=0.28,\qquad y\mid t\sim N(t,0.04^2).
\]
The standard deviation is 4 percentage points, not 4% of the reading. I have chosen the observation for illustration; it is not a reported laboratory result or an extra bulk-assay channel.
This is not a conjugate update. A one-dimensional numerical grid is enough, and makes a useful reference calculation before trying a sampler.
View code
S, y_observed, sigma_y =0.40, 0.28, 0.04def grid_posterior(share_alpha=2., share_beta=2., sigma=0.04, n=4000):if share_alpha <=0or share_beta <=0or sigma <=0or n <100:raiseValueError("Positive prior shapes and sigma, and n >= 100, are required.")# Midpoints avoid evaluating the Beta density at its support boundaries. t = (np.arange(n) +0.5) * S / n log_prior = beta.logpdf(t / S, share_alpha, share_beta) - np.log(S) prior = np.exp(log_prior - logsumexp(log_prior)) log_joint = log_prior + norm.logpdf(y_observed, loc=t, scale=sigma) posterior = np.exp(log_joint - logsumexp(log_joint))return t, prior, posteriordef summarise_grid(t, mass): mean =float(t @ mass)# Interpolate a piecewise-constant density across bin edges. edges = np.linspace(0., S, len(t) +1) lo, hi = np.interp([0.05, 0.95], np.r_[0., np.cumsum(mass)], edges)return mean, float(lo), float(hi)t_grid, prior_mass, posterior_mass = grid_posterior()np.testing.assert_allclose(prior_mass.sum(), 1., atol=atol)np.testing.assert_allclose(posterior_mass.sum(), 1., atol=atol)assay_only_mass = prior_mass * np.ones_like(t_grid)assay_only_mass /= assay_only_mass.sum()np.testing.assert_allclose(assay_only_mass, prior_mass, atol=atol)posterior_mean, ci_low, ci_high = summarise_grid(t_grid, posterior_mass)# Independent quadrature and a finer grid check the numerical result.def unnormalised_density(t):return beta.pdf(t / S, 2., 2.) / S * norm.pdf(y_observed, t, sigma_y)Z = quad(unnormalised_density, 0., S, epsabs=1e-11)[0]quad_mean = quad(lambda t: t * unnormalised_density(t), 0., S, epsabs=1e-11)[0] / Znp.testing.assert_allclose(posterior_mean, quad_mean, atol=1e-7, rtol=0.)fine_t, _, fine_mass = grid_posterior(n=8000)np.testing.assert_allclose(summarise_grid(t_grid, posterior_mass), summarise_grid(fine_t, fine_mass), atol=2e-6, rtol=0.)curve_data = DataFrame({"Anatase (wt%)": np.tile(100* t_grid, 2),"Density per wt%": np.r_[prior_mass, posterior_mass] / (100* S /len(t_grid)),"Distribution": np.repeat(["Conditional prior", "Posterior with phase observation"],len(t_grid)),})posterior_plot = ( ggplot(curve_data, aes("Anatase (wt%)", "Density per wt%", color="Distribution", linetype="Distribution"))+ geom_line(size=1) + theme_minimal()+ theme(figure_size=(7, 4), legend_position="bottom"))display(posterior_plot)print(f"Posterior mean: {100* posterior_mean:.2f} wt%; "f"central 90% credible interval: {100* ci_low:.2f}–{100* ci_high:.2f} wt%")
The new phase-sensitive observation updates the conditional anatase prior. The bulk assay alone leaves this conditional distribution unchanged.
Posterior mean: 27.20 wt%; central 90% credible interval: 21.01–33.25 wt%
The posterior now reflects both the prior split and the new observation. Its 90% credible interval contains 90% of the posterior probability under this model. It is not the LP feasibility range and does not, by itself, establish 90% repeated-sample coverage.
Actually check prior and measurement sensitivity
It is easy to write “results depend on the prior” and move on. Here we can inspect that dependence directly. The first three rows below keep the observation error fixed and change the prior on the anatase share. The last two keep the Beta(2,2) share prior and change the assumed observation error.
All these answers are conditional on quartz and total titanium phase being exactly known. A real analysis should also investigate uncertainty in the assay, the response matrix, the phase list and the observation model. A narrow interval obtained by fixing those quantities is not an estimate of every source of uncertainty.
Posterior prediction is a check, not validation
A posterior draw gives a possible \(t\). To simulate a new response, also draw its measurement error:
Central 90% predictive interval (wt%): [18.01 36.11]
This illustrates the difference between uncertainty in a mineral fraction and uncertainty in a new observation. With one invented, in-sample datum, it is a very limited model check. It does not validate an EMC method against independent mineralogy.
NoteOptional: the same conditional model in NumPyro
The grid is the executed example in this post. The following standalone sampler is optional and is not executed during the article render. It requires NumPyro, JAX and ArviZ in addition to the main environment. Its numerical output has not been validated in this revision.
The important prediction detail is y=None: the fitting call supplies the observation, while the predictive call must release it so a new response can be simulated (NumPyro Contributors, n.d.).
Compare its posterior mean and quantiles with the grid, allowing for Monte Carlo error. Inspect divergences, chain mixing, rank-normalised R-hat and effective sample sizes before interpreting the summaries. PyMC can express the same custom model; a formula interface such as Bambi is more naturally used for a separate regression calibration. A different software package is not a different source of mineralogical evidence.
Taking this back to a real EMC model
The next step is not automatically a bigger algorithm. It is a better-characterised problem: a defensible phase list, consistent analytical bases, appropriate mineral chemistry and paired mineralogical measurements where the chemical system needs them. Evaluate calibration on samples not used to fit or tune it, with domains and deployment conditions in mind. Check mineral errors as well as reconstructed assays; our exact-fit wrong phase list shows why both matter.
Where mineral chemistry is uncertain, fit or propagate it together with abundance. Keep the matching posterior draws together when reconstructing the assay:
\[
\widehat a_s^{(d)}=C_s^{(d)}m_s^{(d)}.
\]
Multiplying separate posterior means generally loses their dependence. Component by component,
Allowing both chemistry and abundance to vary can introduce further confounding. Independent mineral-chemistry measurements and well-justified geological restrictions therefore become more important, not less.
For the toy mixture, the useful result is simple: 60 wt% quartz and 40 wt% total titanium phase, with anatase and rutile unresolved by bulk chemistry. If that total is enough for the decision, EMC has done something useful. If the decision depends on the split, we need a measurement that distinguishes the phases.
That is where I think EMC is strongest: using the assay data we already have, being clear about what they support, and directing the next mineralogical measurement towards something that will actually change the answer.
NoteReproducing this article
From the repository root, run uv sync --frozen, then uv run python scripts/check_emc_article.py. The checker executes the main Python cells in order and deliberately skips the optional NumPyro cell. Use uv run quarto render posts/element-to-mineral-conversion/index.qmd to render the article when Quarto is installed.
The main cells use the repository’s existing NumPy, SciPy, pandas, Plotnine and IPython dependencies. There are no external data reads or network calls. The mixture and phase observation are invented; the predictive draws have a fixed seed. Code is folded to keep the argument readable, and execution errors stop the render rather than appearing as article output.
Working with assay data and limited mineralogy?
Get in touch to discuss whether your conversion method supports the decision you need to make.
Please avoid sending confidential datasets in an initial enquiry.
References
Berry, R.F., Hunt, J.A. and McKnight, S.W. (2011). Estimating mineralogy in bulk samples. First AusIMM International Geometallurgy Conference. Proceedings record and abstract.
Johnson, L.J., Chu, C.H. and Hussey, G.A. (1985). Quantitative clay mineral analysis using simultaneous linear equations. Clays and Clay Minerals, 33(2), 107–117. doi:10.1346/CCMN.1985.0330204.
Lund, C., Lamberg, P. and Lindberg, T. (2013). Practical way to quantify minerals from chemical assays at Malmberget iron ore operations: An important tool for the geometallurgical program. Minerals Engineering, 49, 7–16. doi:10.1016/j.mineng.2013.04.005.
Mena Silva, C., Sørensen, B.E., Aasly, K. and Ellefmo, S.L. (2018). Geometallurgical approach to the element-to-mineral conversion for the Nabbaren nepheline syenite deposit. Minerals, 8(8), 325. doi:10.3390/min8080325.
Parian, M., Lamberg, P., Möckel, R. and Rosenkranz, J. (2015). Analysis of mineral grades for geometallurgy: Combined element-to-mineral conversion and quantitative X-ray diffraction. Minerals Engineering, 82, 25–35. doi:10.1016/j.mineng.2015.04.023.
Rodrigues, A.A., Grose, L., Ailleres, L., Halley, S., Escolme, A., Armit, R., Harandi, M. and Cracknell, M. (2023). Pedras: Modal mineralogy approximations from assay data using Bayesian inference. GeoBerlin 2023, conference abstract. doi:10.48380/2r85-w395.
Whiten, B. (2008). Calculation of mineral composition from chemical assays. Mineral Processing and Extractive Metallurgy Review, 29(2), 83–97. doi:10.1080/08827500701257860.
